extfieldsofstudy:
Computer Science
Mathematics

s2fieldsofstudy:
Computer Science
Mathematics

text:
How to Use Undiscovered Information Inequalities: Direct Applications of the Copy Lemma

We discuss linear programming techniques that help to deduce corollaries of non classic inequalities for Shannon's entropy. We focus on direct applications of the copy lemma. These applications involve implicitly some (known or unknown) non-classic universal inequalities for Shannon's entropy, though we do not derive these inequalities explicitly. To reduce the computational complexity of these problems, we extensively use symmetry considerations. We present two examples of usage of these techniques: we provide a reduced size formal inference of the best known bound for the Ingleton score (originally proven by Dougherty \emph{et al.} with explicitly derived non Shannon type inequalities), and improve the lower bound for the optimal information ratio of the secret sharing scheme for an access structure based on the Vamos matroid.

I. INTRODUCTION
We can associate with an n-tuple of jointly distributed random variables its entropy profile that consists of the 2 n − 1 values of Shannon's entropy of each sub-tuple of the given tuple. We say that a point in R 2 n −1 is entropic if it represents the entropy profile of some distribution. The entropic points satisfy different information inequalities -the constraints that characterize the range of admissible entropies for jointly distributed variables. The most known and well understood information inequalities are so-called Shannon-type inequalities that are defined as linear combinations of several instances of the basic inequality I(A : B|C) ≥ 0, where A, B, C are any subsets (possible empty) of random variables.
In 1998 Z. Zhang and R.W. Yeung discovered the first example of an (unconditional) non Shannon-type information inequality, which was a linear inequality for entropies of a quadruple of random variables that cannot be represented as a combination of basic inequalities [4]. Now many other non Shannon-type information are known (including several infinite families of non Shannon-type information inequality), see, e.g., [6], [11], [17].
The works on non-Shannon type information inequalities are focused on the following fundamental question: can we prove some particular property for Shannon's entropy (usually expressed as an equality or an inequality), assuming some specific constraints (which are also given as equalities and inequalities for Shannon's entropy, including Shannon type and non Shannon type universal information inequalities). Not surprisingly, these arguments typically use the technique of linear programming: at the first stage, to derive new informa-tion inequalities (and to prove that they are indeed new) and then, at the second stage, to apply the derived inequalities to some specific problem (e.g., in secret sharing, network coding, and so on). Many applications of this type involve heavy computations: the number of information inequalities in n random variables grows exponentially with n. It seems that in several problems the progress is stopped since the computational complexity of the corresponding linear programs is too high for the modern computers. So the question arises: can we reduce the computational complexity of the relevant problems of linear programming? In this paper we show that in some cases new progress can be made with a reasonable combination of previously known techniques.
Recently Farràs et al., [26], observed that the two stages of the scheme mentioned above (deriving new information inequalities and applying them) can be merged in a unique problem of linear programming: in a specific application of information inequalities, instead of enumerating explicitly all known non-Shannon information inequalities, we can include in a linear program the constraints from which those inequalities are (or can be) derived. The constraints used in [26] were instances of the Ahlswede-Körner lemma, [7]. As a tool of inference of non Shannon type inequalities, the Ahlswede-Körner lemma is essentially equivalent to the Copy Lemma from [17] "cloning" one single random variable, see [20]. In this paper we use a technique similar to [26], but employing the constraints obtained from the Copy Lemma. This permits us to extend the technique: we apply the Copy Lemma to make "clones" of pairs of jointly distributed random variables (instead of copies of individual random variables used in [26]).
To reduce the dimension of the corresponding problem of linear programming, we use symmetries of the problems under consideration. Our usage of symmetries is similar to the ideas proposed in [24], where the authors suggested to use symmetries to reduce the complexity of computing the outer bounds on network coding capacity. For a study of the general problem of Shannon and non Shannon type inequalities for symmetric entropy points we refer the reader to [25].
To illustrate the power of this technique, we discuss two examples. In the first one we give a formal inference of the best known bound for the Ingleton score (originally proven in [17]), with only four auxiliary random variable and three applications of the Copy Lemma. The obtained computerassistant proof is very fast, and the "formally checkable" version of the proof can be reduced to a combination of 167 equalities and inequalities with rational coefficients. In our second example we improve the known lower bound for the optimal information ratio of the secret sharing schemes associated with the Vámos matroid (for the access structure denoted V 0 ).
Though the techniques used in this paper are rather simple, their combination proves to be fruitful. We believe that this approach can be used in other problems of information theory.

Notation:
We denote [n] = {1, . . . , n}. Given a tuple of jointly distributed random variable (X 1 , . . . , X n ) and a set of indices V = {i 1 , . . . , l k } ⊂ [n] we denote by H(X V ) the Shannon entropy H(X i1 , . . . , X i k ). We fix an arbitrary order on nonempty subsets of the set of indices [n] and assign to each distribution X = (X 1 , . . . , X n ) its entropy profile, i.e., the vector of 2 n −1 entropies H(X V ) (for all non-empty V ⊂ [n]).
We say that a point in the space V n = R 2 n −1 is entropic if it represents the entropy profile of some distribution. Following [4], we denote by Γ * n the set of all entropic points. We denote byΓ * n the topological close Γ * n ; its elements are called almost entropic.
We use the standard abbreviations for linear combinations of coordinates in the entropy profile: We extend this notation to almost entropic points.

II. COPY LEMMA
The following lemma was used (somewhat implicitly) in the very first proof of a non Shannon type inequality in [4], see also [18,Lemma 14.8]. The general version of this lemma appeared in [9] and [17]. In [17] it was presented explicitly, with a complete and detailed proof.
Lemma 1 (Copy Lemma): For every tuples of jointly distributed random variables (X, Y, Z) there exists a distribution (A, B, C, C ′ ) such that • the distribution (A, B, C) coincides with the distribution of (X, Y, Z); the distribution of (A, C ′ ) coincides with the distribution of (X, Z), and In what follows we abuse the notation and identify (A, B, C) with (X, Y, Z). That is, we say that the initial distribution (X, Y, Z) can be extended 1 to the distribution (X, Y, Z, Z ′ ), where (X, Z) and (X, Z ′ ) have isomorphic joint distributions, and I(Z ′ ; Y, Z|X) = 0. Following [17], we say that Z ′ is a Y -copy of Z over X, denoted by Z ′ := Y -copy(Z|X). 1 Formally speaking, such an "extension" may require a subdivision of the initial probabilistic space, so technically (A, B, C) and (X, Y, Z) are defined on different probabilistic spaces.
Note that Lemma 1 remains valid if X, Y, Z are not individual variables but tuples of variables.
F. Matúš and L. Csirmaz proposed several extensions and generalizations of the Copy Lemma (polymatroid convolution, book extension, maximum entropy extension, see, e.g., [12], [22], [23]). However, to the best of our knowledge, the "classic" version of the Copy Lemma is strong enough for all known proofs of non Shannon type inequalities All known proofs of (unconditional) non Shannon type information inequalities can be presented in the following style: we start with a distribution (X 1 , . . . , X n ); using the Copy Lemma we supplement to this distribution several new variable (Y 1 , . . . , Y k ); then we enumerate the basic Shannon's inequalities for the joint distribution (X 1 , . . . , X n , Y 1 , . . . , Y k ) and show that a combination of these inequalities (together with the constraints from the Copy Lemma) gives the desired new inequality for (X 1 , . . . , X n ) (the entropy terms with new variable Y j should be canceled). We refer the reader to [17] for many instances of this argument.
In what follows we use a very similar scheme of proof. The difference is that we superimpose constraints (equalities for the coordinates of the entropy profile of (X 1 , . . . , X n )) that are specific for the given specific problem. That is, in each application we prove a conditional non Shannon type inequality (which is valid only assuming some linear constraints for the entropy quantities). Each of these arguments can be theoretically translated in a more conventional style: we can at first prove some new unconditional non Shannon type inequalities, and then combine them with the given constraints and deduce the desired corollary. However, we do not need to reveal and formulate explicitly these unconditional inequalities hidden behind the proof.

III. SYMMETRIES
In this section we discuss a simple formalism of symmetry considerations that helps to reduce the dimension of a linear programs associated with problems concerning information inequalities (assuming that the problem has a particular symmetric structure). A similar formalism (action of permutations on multidimensional distributions and their entropy profiles, information inequalities for distributions with symmetric entropy profiles) was systematically studied in [25].
Let X (n) = (X 1 , . . . , X n ) be a joint distribution of n random variables. We denote by H(X) its entropy profile. Let π : [n] → [n] be a permutation of indices (an n-permutation). This permutation induces a natural transformation of the distribution X, and therefore a transformation of the entropy profile Example 1: If π is the transposition of indices 1 and 2, then the mapping defined above exchanges the entropy values H(X 1 ) and H(X 2 ), the values H(X 1 , X 3 ) and H(X 2 , X 3 ), etc., and does not change the values H(X 1 , X 2 ), Thus, every permutation of indices π induces a transformation of the set of all entropic points Γ * n . This transformation of the entropic points can be naturally extended to a linear transformation of the entire space V n (defined by a suitable permutation of the coordinate axis). In what follows we denote this transformation by T π . (Note that π is a permutation on n elements, while T π is a transformation of the space of dimension 2 n − 1.) In the usual algebraic terminology, we can say that we defined a representation of the symmetric group (the group of all n-permutations) by linear transformations of the space V n .
We say that a point We denote the set of G-invariant points in V n by Inv(G). In what follows we typically study linear and affine spaces that are G-invariant for some specific group of permutations G. Note that every G-invariant linear space in V n contains the set Inv(G) as a subspace.
The sets of all entropic and all almost entropic points (Γ * n andΓ * n respectively) are obviously π-invariant for every npermutation π. For every group G of n-permutations, the intersections of Γ * n orΓ * n with any G-invariant set gives another G-invariant set. In what follows we discuss several simple examples of this type.
Example 2: Let L ent be the set of points in Γ * 4 such that Since L ent is defined by affine conditions, it can be represented as the intersection of Γ * 4 with an affine subspace L (of co-dimension 5) in the entire space V 4 = R 15 . It is easy to see that this subspace L is invariant with respect to transpositions π 12 and π 34 (π 12 exchanges the indices 1 and 2, and π 12 exchanges the indices 3 and 4). Therefore, this subspace is G-invariant, where G is the subgroup of the symmetric group generated by the transpositions π 12 and π 34 (this subgroup consists of four permutations).
We say that a linear function Similarly, for a subgroup of n-permutations G, we say that a linear function is G-invariant, if it is π-invariant for every π ∈ G.
Example 3: Ingleton's quantity Ing(X 1 , X 2 , X 3 , X 4 ) defined as can be extended to a linear function on the space V 4 = R 15 . This function is invariant with respect to the transposition of indices 1 and 2, to the transposition of indices 3 and 4, and with respect to the group (of size four) generated by these two transpositions.
Claim 1: Let G = {π 1 , . . . , π k } be a subgroup of npermutations, and let x be a point in V n . Then the center of mass of the points is G-invariant. Proof: It is easy to see that the set of points (1) is πinvariant for each π ∈ G. Therefore, the center of mass of this set is also G-invariant.
Lemma 2: Let G = {π 1 , . . . , π k } be a subgroup of npermutations, L be a G-invariant convex set in V n , and F : V n → R be a G-invariant linear function. Then Proof: It is enough to prove that for every x ∈ L there exists a y ∈ L ∩ Inv(G) such that To this end we take the points and define y as the center of mass of these k points. Since L is G-invariant, each point T πi (x) belongs to L. From convexity of L it follows that y also belongs to L, and Claim 1 implies that y belongs to Inv(G). Thus, the constructed point belongs to L ∩ Inv(G).
On the other hand, from G-invariance of F it follows that Since F is linear, we conclude that F at the center of mass of T πi (x) is equal to the same values as the points in (2), and the lemma is proven. Example 4: Assume we are looking for the minimum of Ing(X 1 , X 2 , X 3 , X 4 ) (see Example 3) on the set of almost entropic points (points inΓ * 4 ) that satisfy the linear constraints (see Example 2). We define G as the group of permutations generated by two transpositions π 12 and π 34 that exchange the indices (1, 2) and (3, 4) respectively.
We claim that the required extremal value of Ing(X 1 , X 2 , X 3 , X 4 ) is achieved at some G-invariant point. Indeed, the given constraints define in V 4 a G-invariant affine subspace. The intersection of this subspace with the convex coneΓ * 4 is a G-invariant convex set. On the other hand, Ingleton's quantity Ing(X 1 , X 2 , X 3 , X 4 ) is a G-invariant linear function. Therefore, we can apply Lemma 2 and conclude that desired extremum is achieved at some point in Inv(G).

IV. THE STANDARD BENCHMARK: BOUNDS FOR INGLETON'S QUANTITY
It is widely believed that one of the most interesting questions on the structure ofΓ * 4 is the question on the Ingleton's quantity: what is the worst possible violation of Ingleton's' inequality Ing(A, B, C, D) ≥ 0.
This inequality is true for linearly representable entropic points but for some other entropic points it is violated (see [5]). The question is how far below 0 can go the Ingleton score Ing(A, B, C, D) H(A, B, C, D) (here we use the terminology from [17]). From the Shannon type inequalities it follows only that this score is greater than −1/4. The first non Shannon type inequality from [4] implies that this score is greater than −1/6. The best bound was discovered by Dougherty et al.: Theorem 1 (originally proven in [17] The objective function Ing(A, B, C, D) and the normalization constraint are invariant with respect to the transpositions π A,B and π C,D that exchange the variables A, B and C, D respectively. Therefore, we can apply Lemma 2 and conclude that the required minimal value can be found in the space of points that are π A,B -and π C,D -invariant.
The invariance means that we can restrict ourselves to the points inΓ * 4 such that     Technically, this means that we take the conditions (equalities) from the Copy Lemma (Lemma 1) that determine the entropies of the introduced variables and the properties of conditional independence for each instance of the Copy Lemma. At last, we add to the normalization condition H(A, B, C, D) = 1 and the constraints (4) and (5)  This problem is easily for the standard linear programming solvers (we made experiments with the solvers [27]- [29]). We can extract from the dual solution of this linear program a rational combination of inequalities that gives the desired bound (without floating point arithmetics nor rounding). We present this "formally checkable" proof in Appendix.

V. SECRET SHARING SCHEMES ON THE VÁMOS MATROID
The notion of a secret sharing scheme was introduced by Shamir [2] and Blakley [1]. Nowadays, secret sharing is an important component of many cryptographic protocols. We refer the reader to [15] for an excellent survey of secret sharing and its applications.
The aim of perfect secret sharing is to distribute a secret value among a set of parties in such a way that only the qualified sets of parties can recover the secret value; the nonqualified sets should have no information about the secret. The family of qualified sets of parties is called the access structure of the scheme. This family of sets must be monotone (a superset of a qualified set is also qualified).
The information ratio of a secret sharing schemes is the maximal ratio between the size (the entropy) of a share and the size (the entropy) of the secret. The optimization of this parameter is the central problem of secret sharing. This question was extensively studied for several families of access structures. The main technique for proving lower bounds on complexity of secret sharing schemes is nowadays the linear programming applied to the constraints formulated in terms of classic and non-classic information inequalities, see [19].
In what follows we apply our technique (the Copy Lemma combined with symmetry considerations and linear programming) to one of the classic access structures -the access structure defined on the Vámos matroid (more technically, one of two access structures that can be defined on the Vámos matroid). We skip the general motivation (the Vámos matroid is particularly interesting as a simple example of a non-representable matroid, see [15] for details) and formulate directly the problem of the information ratio for this access structure in terms of Shannon's entropy of the involved random variables.
We define an access structure with parties {1, . . . , 7}. The minimal qualified sets are the 3-sets {1, 2, 3}, {1, 4, 5} and all 4-sets not containing them, with three exceptions {2, 3, 4, 5}, {2, 3, 6, 7}, {4, 5, 4, 7}. We deal with jointly distributed random variables (S 0 , . . . , S 7 ), where S 0 is the secret and S 1 , . . . , S 7 are the shares given to the seven parties involved in the scheme. We require that (i) for every minimal qualified set V listed above . Under these constraints, we are looking for the extremal value Since the set of almost-entropic points is a closed convex cone, we can add the normalization H(S 0 ) = 1 and rewrite the objective function as It is known that the required minimum is not greater than 4/3, [14]. The history of improvements of the lower bound for this problem is presented in the Notice that the objective function, the normalization condition H(S 0 ) = 1, and the linear constraints in (i) and (ii) above are invariant with respect to the permutation that • swap the indices 2 and 3, • swap the indices 4 and 5, • swap the indices 6 and 7, • swap the pairs (2, 3) and (4, 5). Therefore we can apply Lemma 2 and reduce the problem to those almost entropic points that satisfy the symmetries defined above.
These two conditions can be merged in one (more symmetric) constraint Thus, we obtain a linear program with the following constraints: • the conditions (i) and (ii) that define the access structure of the secret sharing scheme, • the normalization H(S 0 ) = 1, • the equalities that follow from the symmetry of the access structure, • the conditions of the two instances of the Copy Lemma, • and the Shannon type inequalities for twelve random variables (S 0 , . . . , S 7 , V ′ , V ′′ , W ′ , W ′′ ) (i.e., for 2 12 − 1 coordinates of their entropy profile). The goal is to optimize the objective function (6).
Remark 2: We do not need to include in the linear program the variables V and W and their entropies -these two variables are used only as a notation in the Copy Lemma. With a help of a computer (we used the solvers [27]- [29]) we find the optimal solution of this linear program: it is equal to 561/491. Thus we establish the following statement.
Theorem 2: The optimal information rate of the secret sharing scheme for the access structure V 0 on the Vámos matroid is not less than 561/491 = 1.142566 . . .

Remark 3:
The standard linear programming solvers use the floating point arithmetic, and we should keep in mind the rounding errors. Our linear program contains only constraints with integer coefficients, so the values of the variables for the optimal solution of this linear program (the primal and the dual) can be taken rational. These rational values can be found by exact computing, without rounding. To compute the exact solution of a linear program we use the rational linear programming solver QSopt_ex, see [10], [27].
The found exact dual solutions (rational linear combinations of constraints) for this linear program consist of more than 10 3 equalities and inequalities. The found rational linear combinations of equalities and inequalities provide a formal mathematical proof of Theorem 2, though such a proof hardly can be called "human-readable." Computing the exact solution of a linear program is rather time consuming. If we only need a lower bound for the optimal solution, it is enough to find a feasible solution of the dual problem that approximates (well enough) the optimal solution. Such an approximation can be found much faster than the exact ration solution, with more conventional floating point LP solvers.

Remark 4:
The linear program constructed in this section would have the same optimal value if we omit the symmetry conditions (note that, unlike the previous section, in the proof of Theorem 2 the Copy Lemma is applied in symmetric settings). However, these symmetry conditions are not useless. In our experiments with exact computations, the symmetry constraints help to find a rational solution of the dual problem with smaller denominators. Also, in several experiments with floating point solvers, the symmetry constraints make the computation (or approximation) of the optimal solution faster. However, this is not a universal rule: in some experiments the symmetry constraints make the computation only slower -it depends on the used linear programming algorithm and on the chosen scaling of the objective function and of the constraints.

VI. CONCLUSION
When we derive balanced information inequalities, the technique of the Copy Lemma with "clones" of individual random variables is equivalent to the technique of the Ahlswede-Körner lemma, see [20]. In this paper we used the Copy Lemma in a stronger form, making "clones" of pairs of correlated variables. In the proofs of Theorem 1 and Theorem 2 it was crucial that we can "duplicate" in one shot a pair of correlated random variables, not only a single random variable The Ahlswede-Körner lemma used in [26] has a clear intuitive meaning: it can be interpreted in terms of extraction of the common information, see [7] and [6]. It would be interesting to find a structural property of the probability distribution that would give an intuitive explanation of the efficiency of the Copy Lemma with two (or more) copied variables.

APPENDIX
In what follows we present a list of 167 equalities and inequalities. An integer linear combination of these inequalities implies Theorem 1.


Large-Scale Outdoor SLAM Based on 2D Lidar

: For autonomous driving, it is important to navigate in an unknown environment. In this paper, we propose a fully automated 2D simultaneous localization and mapping (SLAM) system based on lidar working in large-scale outdoor environments. To improve the accuracy and robustness of the scan matching module, an improved Correlative Scan Matching (CSM) algorithm is proposed. For efﬁcient place recognition, we design an AdaBoost based loop closure detection algorithm which can efﬁciently reject false loop closures. For the SLAM back-end, we propose a light-weight graph optimization algorithm based on incremental smoothing and mapping (iSAM). The performance of our system is veriﬁed on various large-scale datasets including our real-world datasets and the KITTI odometry benchmark. Further comparisons to the state-of-the-art approaches indicate that our system is competitive with established techniques.

Introduction
Simultaneous localization and mapping (SLAM) is one of the key problems for an autonomous mobile robot and has been actively researched for over two decades [1]. As a fundamental module of mobile robot, SLAM has wide application domains ranging from exploration to autonomous driving. While other sensor modalities such as IMU, camera, GPS or RGBD sensors are commonly used, SLAM based on lidar sensor has been very popular, in part because lidar sensor has the advantages of high precision and strong resistance to interference, and is robust to illumination variations. In particular, in this work we focus on large-scale outdoor SLAM based on 2D lidar sensor. 2D lidar sensor enjoys the merits of both high-accuracy range measurement and low-cost compared to the expensive 3D lidar scanners, making it a more suitable choice for cost-sensitive application scenarios.
Typically, a lidar SLAM system consists of a lidar odometry front-end and a back-end that performs the pose graph optimization to obtain a globally optimized map and the full trajectory of the robot. Besides, loop closure detection plays an important role that provides constraints between poses that were captured at the same place but at different times. For the front-end, Iterative Closest Point (ICP) [2], Normalized Distribution Transform (NDT) [3] and Correlative Scan Matcher (CSM) [4] are typical algorithms that were commonly used for 2D lidar scan matching. In particular, ICP and NDT are efficient local scan matching algorithms that depend on an accurate initial guess, but may get stuck in a local optimum if the initial guess is far from the ground-truth. On the contrary, CSM is a global scan matching algorithm that could find the global optimum at a fixed resolution in a probabilistic framework. Various approaches [5,6] have been proposed to improve the performance of CSM for 2D lidar scan matching. Previous works have obtained satisfactory results in relatively small indoor environment, whilst the application in large-scale outdoor environment is still a challenging task.
The primary challenges in designing a 2D lidar SLAM system for large-scale outdoor environment lie in three aspects: Firstly, the front-end should be designed computationally-efficient and robust to initialization error, generating a good prior for loop closure detection and global optimization, especially in low-texture and dynamic environment. Secondly, an efficient loop closure detection algorithm is needed to reduce accumulation drift when loops are detected. This is of particular importance in large-scale settings. However, most existing loop closure detection algorithms are designed for visual SLAM systems [7][8][9][10], and little attention has been given to 2D lidar SLAM systems [11][12][13]. Finally, an accurate and robust optimization back-end is needed to reject false loop closures.
To tackle the above challenges, this paper makes the following three contributions: (1) Firstly, a modified CSM algorithm is proposed to improve the accuracy and robustness of the front-end scan matching algorithm, especially in low-texture and dynamic environment. (2) Secondly, we propose an AdaBoost based loop closure detection algorithm and a false loop closure rejection algorithm which work efficiently to perform place recognition. (3) Thirdly, we propose a light-weight back-end optimization algorithm that works in real-time.
The optimization results are also utilized to eliminate false loop closures.
The remainder of this paper is organized as follows. In Section 2, some related works are introduced. The details of the implementation of our approach are presented in Section 3. Section 4 provides the experimental results and comparisons with state-of-the-art approaches on both publicly available datasets and our own collected datasets. Finally, the conclusions and future works are summarized in Section 5.

Scan Matching Approaches
The aim of scan matching [14] is to estimate the rigid-body transformation between two laser scans. The Iterative Closest Point (ICP) [2] and Iterative Closest Line [15] are the most popular scan matching algorithms which estimate the relative transform by finding the nearest points between the two laser scans. Zhang et al. [16] proposed a low-drift and real-time lidar odometry and mapping (LOAM) algorithm, which selects feature points based on the roughness of the point's local region. The feature points are then matched using either point-to-line or point-to-plane distances. There have been other scan matching approaches, including NDT [3], hough transform-based approaches [17] and histogram-based approaches [18,19]. All these approaches rely on some model assumptions, and may get stuck in a local optimum, especially in large-scale, low-texture and dynamic environments.
Olson et al. [4] proposed the Correlative Scan Matcher (CSM) to search for the best rigid transformation in a probabilistic framework. Compared with other scan matching approaches, CSM is robust to large initialization errors while running in real-time by utilizing the multi-resolution lookup-table rasterization technique. Chong et al. [5] proposed the Extended Correlative Scan Matcher (ECSM) which improves the scan matching accuracy with surface normal information.

Loop Closure Detection
There have been works based on Monte Carlo localization [20] and Rao-Blackwellized filtering [21] to tackle the 2D lidar loop closure detection problem. Recent works on loop closure detection can be roughly classified into three categories: feature based approaches, optimization based approaches and learning based approaches [13].
Bosse et al. [22] proposed to use orientation histogram and projection histogram for loop closure detection. Himstedt et al. [23] augmented the original Fast Laser Interest Point Transform (FLIRT) [24] descriptor with the orientation and distances of landmarks to accelerate place recognition. The approach of cartographer [12] firstly built submaps and then aligned scans to submaps to generate closed loops by a branch-and-bound approach. Li et al. [13] utilized deep learning to perform loop closure detection in indoor environment. The most similar work to our approach is the method proposed in [11], which treat the loop closure detection problem as a classification problem, and use AdaBoost to detect loop closures.

Methodology
As shown in Figure 1, our full system mainly contains three parts: The front-end (Section 3.1), the loop closure detection part (Section 3.2) and the back-end (Section 3.3). The front-end is in charge of key frame selection, data preprocessing, point cloud feature extraction and scan matching. The loop closure detection part is responsible for the detection and validation of the loop closures. The back-end optimizes vehicle poses and global map based on the loop closure detection results and the scan matching results obtained from the front-end. Unlike traditional descriptions, we separate the loop closure part from the front-end to emphasize its importance, especially in mapping large-scale outdoor environment.

Probabilistic Formulation and CSM Overview
Given two lidar scans, CSM is applied to calculate the optimal rigid body transform between the two lidar poses. The scan matching problem can be described in a probabilistic framework. Let the vehicle move from x i−1 to x i , with motion u obtained from the IMU and the wheel encoder. The goal of scan matching is to find the maximum posterior probability p(x i |x i−1 , u, m, z), where m is the environment model and z is the current observation. According to Bayes rule, we have Supposing every lidar observation z j is independent, the observation model can be written as In the context of scan matching, we set the target scan as P, and the source scan as Q. Let T be the rigid body transform that transforms points in Q into P's coordinate. The observation model can be written as Based on Equations (2) and (3), the problem of finding the optimal transformation T * is thus equivalent to search for the maximum log likelihood of the observation model where P raster stands for the grid map rasterized from the target scan P and Q i (i ∈ N) is a single point of scan Q. T corresponds to a 3*3 homogeneous matrix, defined as

The Improved Correlative Scan Matcher
The multi-resolution CSM proposed in [4][5][6] is described in Algorithm 1. This approach could obtain satisfactory results in most structured environment; however, it will encounter difficulties in handling large-scale, low-texture, and dynamic environment, especially for corridors and tunnels. Based on this algorithm, we make the following improvements:

Algorithm 1 The multi-resolution CSM algorithm
Require: Target scan P, Source scan Q, Initial guess T 0 1: Rasterization of lookup table: L high , L low 2: Initialization of scan matching score and transformation: S best (1) Firstly, we improved the rasterization of the lookup table. Instead of rasterizing the single target scan, we utilize historical multi-frame to generate a more robust and more comprehensive rasterized cost table. The multi-frame rasterization algorithm is summarized in Algorithm 2, and the comparison with traditional rasterization approach is shown in Figure 2. To increase the accuracy of the raster table, a weighted sum of historical raster tables can be used, as shown in line 7 of Algorithm 2. A scan that is far away from the current scan will be set a smaller weight. Transform historical scan χ i to coordinate t − 1: 4: for j = 1 to n i do 6: Rasterization of χ i : 7: end if 11: end for 13: end for 14: return L(χ t−1 ) (2) Secondly, the covariance of the transformation T is utilized to recognize scenes that have high matching uncertainty, such as corridors or tunnels. Similar to [4], we calculate the mean and the covariance for T as follows: Suppose T i is a particle corresponding to a single transform along with its weight S i , {T i , S i } is the posterior distribution of T . Let ∑ be the covariance of transformation T , we have Let ∑ x , ∑ y , ∑ z be the diagonal values of ∑ , we define ∑ l and ∑ s as We further define k as the ratio of ∑ s and ∑ l The comparative results of the k-value in different scenes are shown in Figure 3. As shown, scenes that are well constrained in both x and y directions have a value of k nearby 1.0. And scenes like corridors have small values of k, which implies high uncertainty of translation. Figure 4 shows examples taken from the KITTI odometry benchmark [32]. The score distribution is visualized according to the respective translation for a fixed orientation, and the corresponding raw scans are visualized above.   Our improved CSM algorithm utilizes k as the penalty factor ω(k) to increase the odometry accuracy by means of motion model p(x i |x i−1 , u) as where T j − u means the difference between the observation model and the motion model, and T * corresponds to T j with the maximum score S j . We set ω(k) = λ(k − 1), where λ > 0.

Loop Closure Detection and Validation
Our loop closure detection and validation method is mainly composed of three parts: Detecting loop closure candidates based on AdaBoost, hierarchically selecting positive loops with high confidence and rejecting false loops according to the residual error from the back-end optimizer. It should be noted that, our method of detecting loop closure candidates is similar to [11]. However, we improve the accuracy of classification by introducing more robust features and decrease the false alarm by several validation steps.

Point Cloud Feature Extraction
Similar to [11], we extract rotation invariant features as follows: (1) Area features: Area, close area; (2) Distance features: Distance, close distance, far distance, average range, centroid, maximum range, mean deviation, standard deviation of range; (3) Shape features: Circularity radius, circularity residual, curvature mean, curvature standard deviation, regularity; (4) Cluster features: Number of groups, mean group size; and (5) Others: Size and mean of angular difference. For detailed definition and formulation of these features, please refer to [11,33]. Based on [11,33], we have improved the feature extraction in three aspects.
(1) Firstly, for the mean and the standard deviation of the curvatures, the original method in [11] only takes three consecutive points into account to calculate the curvature of the middle point, which is sensitive to noise. Let p i be a point in the lidar cloud P k , and let S i be the neighboring consecutive point sets of p i . We choose three points on each side of p i , and the smoothness is evaluated as The curvature mean c and the standard deviation σ(c) are defined as (2) Secondly, for the number of groups and the mean of group size, rather than using a constant maximum distance to separate consecutive points [11], we adopt an adaptive breakpoint detector as introduced in [34]. The adaptive distance threshold of point p i is defined as where r i is the distance from p i to the origin of the coordinate, ∆ϕ corresponds to the angle interval of the range scan, λ is an empirical value determining the radius of the threshold circle and σ corresponds to the measurement accuracy of the laser scanner. A breakpoint is detected when the distance between consecutive points exceeds D i Th . An example is shown in Figure 5, where the black points represent breakpoints and the red point sets between two consecutive breakpoints belong to one cluster. Besides the above two features, we introduce another five-dimensional feature based on the size of the cluster. We discretize the size of the cluster into 5 levels, and calculate the ratio of each level L = [l 1 , l 2 , l 3 , l 4 , l 5 ], where l i is normalized. We termed this feature as Cluster Distribution Histogram (CDH). The size of each level is set empirically. This five-dimensional feature together with the first two features constitutes the distribution of the point cloud clusters. Three examples of CDH from our real scan data is shown in Figure 6. (3) Thirdly, we introduce three more features based on scan matching result, including the normalized score of the maximum scan matching alignment, the k-value computed by Equation (7) and the missing rate of query scan points. In particular, the k-value reflects the characteristic of the scene, and the missing rate corresponds to the rate of the query scan points which are projected into the black region of the rasterized cost table as shown in Figure 2.

Classification Based on AdaBoost
All those m rotation invariant features are computed for each scan, defined as where N corresponds to the total number of scans. Given two laser scans k and l, the set of features for the classifier is defined as For the training data, a set of n labeled scan pairs are provided where y i ∈ {1, −1} is the label of a sample input indicating matching and nonmatching scan pairs. AdaBoost is applied as the classifier constituted with a series of weak classifiers. The AdaBoost algorithm is described in Algorithm 3. The algorithm firstly trains a base classifier from the initial training sets, and then adjusts the distribution of the training samples according to the classification error, so that the previously misclassified samples receive more attention in the following step. And then the next weak classifier is trained based on the updated sample distribution. The above two steps are repeated until the maximum iteration T is reached. The final classifier is a weighted combination of T weak classifiers, and these weak classifiers with smaller classification errors have larger weights.

Loop Closure Validation
It is worth noting that one single wrong loop closure can cause severe damage to the whole pose graph. It is necessary to perform validation to eliminate the false loop closures. Our system adopts three strategies to reject false loops: (1) Firstly, for every positive loop closure detected by the AdaBoost classifier, we calculate the best score of the query scan projected to the coarse raster table, defined as S best L . The candidates detected from the classifier with S best L < S Th L are rejected. (2) Secondly, we calculate the scan matching results at high resolution, including the best score S best H , the k computed by Equation (8) and the missing rate of query scan points defined as p miss . Three thresholds are used to select positive loop closures with high confidence.
(3) Thirdly, we perform consistency validation by inverse scan alignment which we named as Direct-Inverse Matching Consistency (DIMC). Unlike local matching algorithms such as ICP and ICL, CSM may produce a large deviation between the direct matching and the inverse matching when two scans belong to a false loop closure. An example from our data is shown in Figure 7, where the blue points belong to the target scan and the red points belong to the source scan. Due to the significant changes in viewpoint, the scan matching result is unreliable, which can be verified through DIMC.

Algorithm 3 Classifier based on AdaBoost
Require: Training set D = {(x 1 , y 1 ), (x 2 , y 2 ), . . . , (x n , y n )}, Maximum iteration T 1: Initialize the sample weight distribution: ω 1 (x) = 1 n ; 2: for t = 1 to T do 3: Train a weak classifier h t based on ω t : Compute classification error of h t : is the authenticity discriminant function. 5: if ε t ≥ 0.5 then 6: continue; 7: end if 8: Compute the weight of the weak classifier h t : Update the sample weight distribution: where η t is the normalization factor. 10: end for 11:

Back-End Optimization
The back-end performs graph optimization based on the results generated from the front-end and the loop closure detection, resulting in a consistent trajectory and an optimized map.
The problem can be represented as the Bayesian belief network [25]. The joint probability of the variables is modeled as where X is the set of nodes and U stands for the set of edges, P(x i+1 |x i , u i ) is the motion model constructed from adjacent frame constraints given by the front-end, P(x j |x i , u ij ) corresponds to the measurement model generated from loop closure detections. The state transformation model corresponding to P(x i+1 |x i , u i ) is defined as where f i,i+1 (·) is the state transformation function, w i,i+1 is the Gaussian noise with zero-mean and covariance ∑ i . Similarly, the state transformation model corresponding to P(x j |x i , u ij ) is defined as and the covariance of w i,j is Λ ij . Combining ( (15), (16), (17)), the maximum-a-posteriori estimation can be formed into a nonlinear optimization problem where e 2 ∑ = e T ∑ −1 e. As described in [25], the above nonlinear optimization problem can be transformed into the standard least squares problem based on Gauss-Newton or Levenberg-Marquardt algorithm where X ∈ R n contains all pose variables, the sparse measurement matrix A ∈ R m×n contains all edges, and b ∈ R m is the right hand side (RHS) vector corresponding to the difference between the estimated state and the measurement vector. Different from the work in [25], we improved the back-end mainly in two aspects.
(1) Firstly, the covariance matrix computed in Equation (6) can be embedded into the above optimization framework as an adaptive information matrix. This will ensure that edges with strong confidence are not prone to be adjusted during the optimization procedure, whilst edges with low confidence are more likely to be tuned.
(2) Secondly, our system rejects false loops according to the residual error from the back-end optimizer. The residual error is computed between trajectories before and after the optimization. A large difference indicates a high likelihood of false loops.

Experiments and Results
During our experiments, the algorithms run on a laptop computer with 2.4 GHz Intel Core i7 CPU and 8 GB of RAM. In order to achieve real-time performance, multiple CPU threads are utilized in parallel. We perform experiments both on the data collected by our own platform and also on KITTI benchmark. Our platform is a golf cart equipped with a Velodyne VLP-16, a wheel encoder and a low cost IMU, as shown in Figure 8. Although VLP-16 produces 16 scanning lines, we only utilize one single horizontal scanning line to mimic a 2D lidar. The same method is also applied on the KITTI datasets.

The Results of the Improved Correlative Scan Matcher
For the first experiment, we perform the simple laser scan alignment without loop closure detection or back-end optimization to validate the improvement of our front-end method. The comparative results of the original CSM and our improved CSM are shown in Figure 9. The translation errors are evaluated following the KITTI's protocol with lengths (100 m, 200 m, . . . , 800 m) in 2D coordinates [32].
As shown in Figure 9, the trajectories of our improved CSM are closer to the ground truth than those of the original CSM in both the 0th sequence and the 5th sequence. While the average running time of our improved CSM increases by less than 0.01 s.

The Results of Loop Closure Detection and Validation
As described in Section 3.2, the loop closure classifier based on AdaBoost is trained from our dataset campus #1. The dataset is collected around the campus covering an area of 1750 m × 910 m while traveling around 15.99 kilometers with 22,935 scans. We choose 2000 matching and 2000 nonmatching laser pairs which are pre-labeled to train the classifier by 10-fold cross validation.
We choose three other approaches for comparison: Method 1 [11] is the most similar one to our method, and the main difference is that we utilize more features as described in Section 3.2; method 2 is based on Support Vector Machine (SVM) using the same features as our method; method 3 is based on deep learning algorithm similar to [13]. Figure 10 shows the precision-recall curves for all these three approaches together with our proposed approach.
Compared with method 1, our algorithm performs better on the campus #1 dataset. A threshold of K = 0.3 gives a detection rate of 92.6% and a false alarm rate of 2.7%. Although the false alarm rate of our method is a bit higher than method 1, it will be further reduced through validation (DIMC) and false loop closure rejection based on the back-end optimization. The subsequent experiments will show that our system can achieve nearly 0% false alarm for different types of datasets. Method 3 achieves the worst performance on our dataset. The main reason is that the dataset tested in [13] is a much smaller dataset collected in structured indoor environment whilst our campus #1 dataset is built in large-scale outdoor environment. For time complexity, our method outperforms method 1 both in feature extraction and classification, as shown in Table 1. During all the experiments including fifteen datasets from campus #1, one dataset from the Zhongdian Science and Technology Park and eleven datasets from KITTI benchmark, only two false loop closures are rejected according to the residual error from the back-end optimizer. Our loop closure detection and validation algorithm performs consistently well. Figure 11 shows the mapping of campus #1 area, where the map (shown in red dots) is projected onto a satellite image. The map shows a good consistency with the satellite image. Besides campus #1, we also test our 2D-SLAM system on several other datasets including a dataset from the Zhongdian Science and Technology Park as shown in Figure 12. Results indicate that our system works consistently well in large-scale outdoor environment. Figure 13 shows a snapshot of our mapping process. The red points are obstacles mapped online; the blue points represent the current laser scan; the green points represent the vehicle trajectory and the dark yellow lines are the loop closures. The current vehicle pose is shown as the yellow arrow.   We qualitatively compare our method with a state-of-the-art 2D SLAM algorithm: Cartographer [12]. Figure 14 shows the trajectories generated by our method and cartographer. As can be seen from Figure 14, our method can compete with cartographer in terms of accuracy. For computational efficiency, our method built the whole map in 27.54 minutes compared with 31.09 min for cartographer. Detailed comparisons of these two approaches are shown in Figure 15, where scenes #1, #2 and #3 are labeled in Figure 14. For scene #1, cartographer generates a wrong loop closure which results in blurry lines at place A. As for our map rendering step applies dynamic objects filtering borrowed from Wang et al. [35], the unreliable points are filtered effectively. The different results can be seen from B and C in scene #1, A and B in scene #2 as well as C in scene #3. Place C in scene #2 along with A and B in scene #3 indicate that our map contains richer details with less noise than cartographer.

Test with KITTI Datasets
Our system has also been tested on the large-scale and publicly available KITTI datasets, and the experimental results are shown in Figure 16. To mimic the 2D laser setting, we extract one single laser scan of the Velodyne HDL-64 lidar. Due to space issue, we present results of two typical sequences (Sequences 13 and 15). We choose LOAM, ICP and NDT as the comparison approaches. ICP and NDT are widely used scan matching approaches in lidar SLAM system. LOAM is ranked #1 among all lidar based methods on KITTI dataset.
It is worth mentioning that all these comparison approaches are based on the whole 3D lidar data, whilst our approach only utilizes one single scanning line of the 3D lidar. As shown in Figure 16, the experimental results indicate that our system achieves comparable performance compared to other benchmark approaches. It is worth mentioning that sequence 13 is prone to false loop closures around the viaduct due to the limitation of 2D framework. However, our system can detect and eliminate these wrong loops successfully. For computational efficiency, our system achieves a superior runtime of 0.06 s per scan compared with 0.1 s of LOAM.

Concluding Remarks
In this paper, a robust 2D lidar SLAM system for large-scale outdoor environment is proposed. To improve the accuracy and robustness of scan matching, an extension to the CSM is proposed. For efficient place recognition, an AdaBoost based loop closure detection is proposed. We also propose a light-weight back-end optimization algorithm that works in real-time. The optimization results are also utilized to eliminate false loop closures. Experiments both on our own collected dataset and the publicly available KITTI dataset confirmed the effectiveness of our approach. For the future work, we would like to apply our method to 3D mapping and localization. Moreover, the implementation of our map for localization on different platforms will be verified in our further work.


text:
Joint Task Assignment and Resource Allocation for D2D-Enabled Mobile-Edge Computing

With the proliferation of computation-extensive and latency-critical applications in the 5G and beyond networks, mobile-edge computing (MEC) or fog computing, which provides cloud-like computation and/or storage capabilities at the network edge, is envisioned to reduce computation latency as well as to conserve energy for wireless devices (WDs). This paper studies a novel device-to-device (D2D)-enabled multi-helper MEC system, in which a local user solicits its nearby WDs serving as helpers for cooperative computation. We assume a time division multiple access (TDMA) transmission protocol, under which the local user offloads the tasks to multiple helpers and downloads the results from them over orthogonal pre-scheduled time slots. Under this setup, we minimize the computation latency by optimizing the local user's task assignment jointly with the time and rate for task offloading and results downloading, as well as the computation frequency for task execution, subject to individual energy and computation capacity constraints at the local user and the helpers. However, the formulated problem is a mixed-integer non-linear program (MINLP) that is difficult to solve. To tackle this challenge, we propose an efficient algorithm by first relaxing the original problem into a convex one, and then constructing a suboptimal task assignment solution based on the obtained optimal one. Next, we consider a benchmark scheme that endows the WDs with their maximum computation capacities. To further reduce the implementation complexity, we also develop a heuristic scheme based on the greedy task assignment. Finally, numerical results validate the effectiveness of our proposed algorithm, as compared against the heuristic scheme and other benchmark ones without either joint optimization of radio and computation resources or task assignment design.

I. INTRODUCTION
It is envisioned that by the year of 2020, around 50 billions of interconnected Internet of things (IoT) devices will surge in wireless networks, featuring new applications such as video stream analysis, augmented reality, and autonomous driving. The unprecedented growth of these latencycritical services requires extensive real-time computation, which, however, is hardly affordable by conventional mobile-cloud computing (MCC) systems that usually deploy cloud servers far away from end users [2]. Compared with MCC, mobile-edge computing (MEC) endows cloudcomputing capabilities within the radio access network (RAN) such that the users can offload computation tasks to edge servers in their proximity for remote execution and then collect the results from them with enhanced energy efficiency and reduced latency (see [3] and the references therein). Meanwhile, in industry, technical specifications and standard regulations are also being developed by, e.g., the European Telecommunications Standard Institute (ETSI) [4], Cisco [5], and the 3rd Generation Partnership Project (3GPP) [6].
Various efforts have been dedicated to addressing technical challenges against different computation offloading models. In general, computation task models fall into two major categories, namely the partial offloading and the binary offloading. Tasks that cannot be partitioned for execution belong to the former category, while the latter category supports fine-grained task partitions.
Under a single-user single-server setup, [7] jointly optimized the task offloading scheduling and the transmitting power allocation to minimize the weighted sum of the execution delay and the system energy consumption. The optimal communication strategy as well as computational load distribution was obtained in [8] in terms of the trade-offs of energy consumption versus latency for a multi-antenna single-user single-server system. In a single-user single-helper singleserver system, [9] proposed a novel user cooperation in both computation and communication to improve the energy efficiency for latency-constrained computation. In a multiple-input multipleoutput (MIMO) multicell system, [10] jointly optimized the precoding matrices of multiple wireless devices (WDs) and the CPU frequency assigned to each device with fixed binary offloading decisions, in order to minimize the overall users' energy consumption. A multi-user MCC system with a computing access point (CAP), which can serve as both a network gateway connecting to the cloud and an edge cloudlet, was studied in [11] to find the binary offloading decisions. Joint optimization of (partial) offloaded data length and offloading time/subcarrier allocation was studied in a multi-user single-server MEC system based on time-division multiple access (TDMA) and orthogonal frequency-division multiple access (OFDMA), respectively, in [13]. Furthermore, both binary and partial offloading were considered in a multi-user singleserver MEC system exploiting multi-antenna non-orthogonal multiple access (NOMA) in [14].
References [15] and [16] leveraged the inherent collaborative properties of augmented reality (AR) applications across multiple WDs to minimize the users' total energy expenditure.
In the above works, the edge servers are mostly assumed to be one integrated server. However, considering multi-user multi-server systems where more than one CAPs are distributed over the network, it becomes non-trivial to answer the fundamental questions such as how to distribute the tasks among multiple servers, and how to schedule multiple tasks on one single server [19][20][21][22][23]. Computation resource sharing among WDs with intermittent connectivity was considered as early as in [19], in which a greedy task dissemination algorithm was developed to minimize task completion time. A polynomial-time task assignment scheme for tasks with inter-dependency was developed in [21] to achieve guaranteed latency-energy trade-offs. However, this line of works often assumed the communication conditions (e.g., transmission rate and multiple access schemes) and/or computation capacities (e.g., execution rate) to be fixed or estimable by some service profiling, but ignored potential performance improvement brought by dynamic management over such resources (e.g., transmitting power, bandwidth, and computation frequency).
In this paper, we study a device-to-device (D2D)-enabled multi-helper MEC system, in which a local user offloads a number of independent computation tasks to multiple nearby WDs serving as helpers, such as smart wearable devices, cellphones, tablets, and laptops, via direct D2D links. The motivation for us to study efficient task assignment and wireless resource allocation algorithms to facilitate D2D-enabled MEC are two-fold. First, as WDs constantly improve their capabilities (e.g., battery capacity, processing speed, and spectral efficiency), heterogeneity of radio and computation resources among WDs can be exploited to support various demanding services while achieving mutual benefits [20]. Second, with the proliferation of WDs, some of them may be prohibited from directly accessing to BSs. In these cases, they can entrust some virtual central controllers managed by network operators to collect their task features, and assist in D2D-enabled MEC by pooling and sharing the resources among each other. Assuming that the tasks cannot be further partitioned, we consider a TDMA communication protocol, under which the local user offloads tasks to different helpers and downloads computation results from them over orthogonal pre-scheduled time slots. We aim for minimizing the overall latency subject to individual energy and computation capacity constraints at the local user and the helpers.
A task offloading framework, D2D fogging, was proposed in [24], where WDs could share the computation and communication resources among each other via the assistance of network operators, and dynamic task offloading decisions were made to minimize the time-averaged total energy consumption. From the perspective of system model, we employ different communication protocol from that in [24]. We adopt a three-phase TDMA protocol, under which task assignment becomes a generally NP-hard problem because the tasks offloaded to different helpers (scheduled in different TDMA slots) are in couple with each other. By contrast, under OFDMA, all tasks can be executed independent of each other subject to a common deadline constraint as shown in [24]. Furthermore, as in this paper, each WD is assumed to be assigned with more than one task, the efficient matching-based algorithm that forms the building block of the online task assignment scheme in [24] cannot be applied any more. It is also worth noting that the major difference between this paper and the earlier conference version [1] is that instead of fixing the processing capacities, we consider controllable computation frequencies by exploiting dynamic voltage and frequency scaling (DVFS) [25] to achieve improved overall latency. To our best knowledge, this paper is among the earliest works investigating TDMA-based joint binary task offloading and wireless resource allocation for multiple tasks in a single-user multi-helper MEC system.
The contributions of our paper are summarized as follows. First, we transform the computation latency minimization problem with complex objective function into an equivalent one by investigating the optimal structure of the solution. Next, we jointly optimize the tasks assignment, the task offloading time/rate, the local and remote task execution time/computation frequency, and the results downloading time/rate, subject to individual energy and computation frequency constraints at the local user and the helpers. However, since the formulated problem is a mixedinteger non-linear program (MINLP) that is difficult to solve in general, we propose an efficient algorithm to obtain a high-quality sub-optimal solution by relaxing the binary task assignment variables into continuous ones, and then constructing a suboptimal task assignment solution based on the optimal one to the relaxed problem. Furthermore, to reduce the implementation complexity, we also provide fixed-frequency task assignment and wireless resource allocation as a benchmark, and design a greedy task assignment based joint optimization algorithm. Finally, we evaluate the performance of the proposed convex-relaxation-based algorithm as compared against the heuristic one and other benchmark schemes without joint optimization of radio and computation resources or without task assignment design.
The remainder of this paper is organized as follows. The system model is presented in Section II. The joint task assignment and wireless resource allocation problem is formulated in Section III. The convex-relaxation-based joint task assignment and wireless resource allocation algorithm is proposed in Section IV, while two low-complexity benchmark schemes are proposed in Section V. Numerical results are provided in Section VI, with concluding remarks drawn in Section VII.
Notation-We use upper-case boldface letters for matrices and lower-case boldface ones for vectors. "Independent and identically distributed" is simplified as i.i.d., and means "denoted by". A circularly symmetric complex Gaussian (CSCG) distributed random variable (RV) y with mean u and variance σ 2 is denoted by y ∼ CN (u, σ 2 ). A continuous RV z uniformly distributed over [a, b] is denoted by z ∼ U[a, b]. R M ×N and R N stand for the sets of real matrices of dimension M × N and real vectors of dimension N, respectively. The cardinality of a set is represented by | · |. In addition, P(N) means an N-degree polynomial.

II. SYSTEM MODEL
We consider a multi-user cooperative MEC system that consists of one local user, and K nearby helpers denoted by the set K = {1, . . . , K}, all equipped with single antenna. For convenience, we define the local user as the (K + 1)-th WD. Suppose that the local user has L independent tasks 1 to be executed, denoted by the set L = {1, . . . , L}, and the input/output data length of each task is denoted by T l /R l in bits. In the considered MEC system, each task can be either computed locally, or offloaded to one of the K helpers for remote execution. Let Π ∈ R L×(K+1) 1 In this paper, we do not consider interdependency among tasks enabling data transmission from one helper to another as in [19,21], since even under this simple task model, it becomes clear later that task assignment among multiple D2D helpers over pre-scheduled TDMA slots has already been very demanding to solve. denote the task assignment matrix, whose (l, k)-th entry, denoted by π(l, k) ∈ {0, 1}, l ∈ L, k ∈ K ∪ {K + 1}, is given by 1, if the lth task is assigned to the kth WD, 0, otherwise.
Also, define L (k) = {l ∈ L : π(l, k) = 1} as the set of tasks that are assigned to WD k, k ∈ K ∪ {K + 1}. It is worthy of noting that we assume |L (k) | ≥ 1, k ∈ K ∪ {K + 1}. That is, each WD including the local user should be assigned with at least one task, i.e., L ≥ K + 1 2 .
Define by C l in cycles the number of CPU cycles required for computing the lth task, l ∈ L [20,22]. Also, denote the CPU frequency in cycles per second (Hz) at the kth WD as f k ,

A. Local Computing
The tasks in the set L (K+1) are executed locally with the local computation frequency in cycles per second given as [17] 2 In practice, when L ≤ K, a group of K ′ helpers are required to be selected a priori such that L ≥ K ′ + 1. However, detailed design regarding such selection mechanism is out of the discussion of this paper, and is left as our future work.
where t c 0 denotes the associated local computation time, and f 0 is subject to the maximum frequency constraint, i.e., f 0 ≤ f max 0 . The corresponding computation energy consumed by the local user is given by [25] where κ 0 is a constant denoting the effective capacitance coefficient that is decided by the chip architecture of the local user. Replacing f 0 in (2) with (1), E c 0 can thus be expressed in terms of t c 0 as follows:

B. Remote Computing at Helpers
The tasks assigned in L (k) is offloaded to the kth helper, k ∈ K, for remote execution. In this paper, we consider a three-phase TDMA communication protocol. As shown in Fig. 1, the local user first offloads the tasks in the set L (k) to the kth helper, k ∈ K, in a pre-scheduled order 3 via TDMA in the task offloading phase. Then the helpers execute their assigned computation tasks in the task execution phase. At last, in the results downloading phase, the helpers send the results back to the local user in the same order as in the task offloading phase via TDMA. Note that at each TDMA time slot during the task offloading phase, the local user only offloads tasks to one helper. Similarly, during the results downloading phase, only one helper can transmit over each time slot. In the following, we introduce the three-phase protocol in detail.

1) Phase I: Task
Offloading First, the tasks are offloaded to the helpers via TDMA. For simplicity, in this paper we assume that the local user offloads the tasks to the helpers with a fixed order of 1, 2, . . . , K as in Fig. 1. In other words, the local user offloads tasks L (1) to the 1st helper, then L (2) to the 2nd helper, until L (K) to the Kth helper.
Let h k denote the channel power gain from the local user to the kth helper for offloading, k ∈ K. The achievable offloading rate (in bits per second) at the kth helper is given by 3 Since frequent change of TDMA scheduling policy incurs large amount of signalling overhead, we assume a fixed-order TDMA protocol in this paper, which is practically reasonable, and also commonly adopted in the literature, e.g., [13,17].
where B in Hz denotes the available transmission bandwidth, p of f k is the transmitting power for offloading tasks to the kth helper, and σ 2 k is the power of additive white Gaussian noise (AWGN) at the kth helper. Then, the time spent in offloading tasks to the kth helper is given by According to (4) and (5), whereh k = h k /σ 2 k is the normalized channel power gain, and f (x) is a function defined as f (x) 2 x B − 1. The total energy consumed by the local user for offloading all the tasks in L is thus expressed as 2) Phase II: Task Execution After receiving the assigned tasks L (k) , k ∈ K, the kth helper proceeds with the computation frequency given by where t c k 's is the remote computation time spent by the kth helper. Similarly, helper k's remote computing frequency given by (8) is also constrained by its maximum frequency, i.e.,f k ≤ f max k .
In addition, its computation energy is expressed as where κ k is the corresponding capacitance constant of the kth helper.

3) Phase III: Results Downloading
After computing all the assigned tasks, the helpers begin transmitting the computation results back to the local user via TDMA. Similar to the task offloading phase, we assume that the helpers transmit their respective results in the fixed order of 1, . . . , K. Let g k denote the channel power gain from helper k to the local user for downloading. The achievable downloading rate from the kth helper is then given by where p dl k denotes the transmitting power of the kth helper, and σ 2 0 denotes the power of AWGN at the local user. The corresponding downloading time is thus given by Combining (10) and (11), the transmitting power of the kth helper is expressed as whereḡ k = g k /σ 2 0 denotes the normalized channel power gain from the kth helper to the local user. The communication energy consumed by the kth helper is thus given by

C. Total Latency
Since TDMA is used in both Phase I and Phase III, each helper has to wait until it is scheduled.
Specifically, the first scheduled helper, i.e., helper 1, can transmit its task results to the local user only when the following two conditions are satisfied: first, its computation has been completed; and second, task offloading from the local user to all of the K helpers are completed such that the wireless channels begin available for data downloading. As a result, helper 1 starts transmitting its results after a period of waiting time given by where t c 1 is the task execution time at helper 1. Moreover, for each of the other K − 1 helpers, it can transmit the results to the local user only when: first, its computation has been completed; second, the (k − 1)th helper scheduled preceding to it has finished transmitting. Consequently, denoting the waiting time for helper k (k ≥ 2) to transmit the results as I k , I k is expressed as Accordingly, the completion time for all the results to finish downloading is expressed as To sum up, taking local computing into account as well, the total latency for executing all of the L tasks is given by

III. PROBLEM FORMULATION
In this paper, we aim at minimizing the total latency for local/remote computing of all the tasks by optimizing the task assignment strategy (π(l, k)'s), the task offloading time (t of f k 's), the task execution time (t c k 's), and the results downloading time (t dl k 's), subject to the individual energy and computation frequency constraints at the local user as well as the K helpers. Specifically, we are interested in the following problem: In the above problem, the objective function T total is given by (17). The constraints given by (18a) and (18b) state that the total energy consumption of computation and transmission for the local user and the kth helper cannot exceed E 0 and E k 's, respectively. In (18a), E c 0 and E of f 0 are replaced with (2) and (7), respectively, while (18b) is obtained by substituting (9) and (13) for E c k 's and E dl k 's, respectively. (18c) and (18d) guarantee that the computation frequencies of the local users (c.f. (1)) and the helpers (c.f. (8)) stay below their respective limits. (18e) guarantees that each task must be and only assigned to one WD; and (18f) ensures that each of the local user and the helpers is assigned with at least one task. Finally, (18g) imposes the binary offloading constraints.

A. Problem Reformulation
Note that T total (c.f. (17)) is a complicated function involving accumulative max(·) mainly due to the recursive expression of I k (c.f. (15)). Hence, to obtain an explicit objective function in terms of the optimization variables, we need to simplify T total exploiting the following proposition.
Proposition 3.1: Problem (P0) can be recast into an equivalent problem as follows.
Proof: A brief idea of the proof is given as follows. To remove max(·) in I k 's, we first need to narrow down from different cases leveraging the property of the optimal solution. Then, based on the simplified case, we recursively derive I k 's for k ≥ 2. Finally, we arrive at a clear objective function of (P1-Eqv) subject to all the optimality conditions given by (19a)-(19c). Please refer to Appendix A for the proof in detail.

B. Suboptimal Design
The transformed problem (P0-Eqv) is seen as an MINLP due to the integer constraints given by (18g), and is thus in general NP-hard. Although the optimal solution to (P0-Eqv) can be obtained by exhaustive search, it is computationally too expensive (approx. O((K + 1) L ) times of search) to implement in practice. Therefore, we solicit two approaches for suboptimal solution to (P0-Eqv) in the following sections. The first approach is to relax the binary variables into continuous ones while the second approach aims for decoupling the task assignment and wireless resource allocation.
For the first approach, first, we relax (18g) into continuous constraints expressed as Therefore, the relaxed problem is expressed as: (20).
It is worthy of noting that, since E c 0 (c.f. (3)) and E of f 0 (c.f. (7)) are obtained by convex operations on perspective of convex functions ( l π(l, K + 1)C l ) 3 and f ( l π(l, k)T l )'s with respect to (w.r.t.) the variables t c 0 and t of f k 's, respectively, they are also convex functions. So are E c k and E dl k , ∀k ∈ K. Therefore, (P1) is a convex problem. Next, we need to round the continuous π(l, k)'s into binary one such that (18e) and (18f) are satisfied. The details of the proposed joint task assignment and wireless resource allocation scheme will be discussed in Section IV. In addition, we also provide a brief discussion regarding one special case of this approach in Section V-A, in which computation frequencies of all the WDs are fixed to be their maximum, thus serving as a benchmark scheme without computation allocation.
For the second approach, first, it is easy to verify that given Π fixed, (P0-Eqv) reduces to be a convex problem shown as below: Then, we decouple the design of task assignment and wireless resource allocation by employing a greedy task assignment based heuristic algorithm that will be elaborated in Section V-B.

IV. JOINT TASK ASSIGNMENT AND WIRELESS RESOURCE ALLOCATION
The main thrust of the proposed scheme in this section is to relax the binary task-assignment variables into continuous ones, and to solve the relaxed convex problem in semi-closed forms, which are then followed by attaining suboptimal task assignment based on the optimal solution to the relaxed problem.
It is seen that problem (P1) is convex, and can thus be efficiently solved by some off-theshelf convex optimization tools such as CVX [26]. To gain more insights into the optimal rate and computation frequency allocation, in this section, we propose to solve (P1) leveraging the technique of Lagrangian dual decomposition. The (partial) Lagrangian of (P1) is expressed as where η, β 0 , λ 0 , and ζ 0 denote the dual variables associated with the constraints (19a), (19c), (18a), and (18c), respectively; λ = (λ 1 , . . . , λ k ) T represent the dual variables associated with the total energy constraints (18b) each for one helper; β = (β 2 , . . . , β K ) T are the dual variables for the constraints given by (19b); and the multipliers ζ = (ζ 1 , . . . , ζ K ) T are assigned to the constraints given by (18d). After some manipulations, (21) can be equivalently expressed as and with A k , B k , D k , ∀k ∈ K given by and respectively. The dual function corresponding to (22) can be expressed as As a result, the dual problem of (P1) is formulated as

A. Dual-Optimal Solution to (P1)
In this subsection, we aim for solving problem (P1-dual). To facilitate solving the optimum {t of f k , t dl k , t c k } and t c 0 to (27) providing that Π =Π and a set of dual variables are given, we decompose the above problem into K + 1 subproblems including K for ∀k ∈ K and one for k = K + 1 as follows.
Since these problems are independent of each other, they can be solved in parallel each for one k, k ∈ K ∪ {K + 1}.
is the principal branch of Lambert W function defined as the inverse function of xe x = y [27]. Then, in accordance with the optimal solution to the above subproblems, the optimal time and power together with the optimal task assignment to (27) are shown in the following proposition.
In addition,π(l, k)'s shown in (29) and (30) denote the optimal solution to the following linear programming (LP) problem: where φ l,k , ∀k ∈ K, ∀l ∈ L, is given by and φ l,K+1 , ∀l ∈ L, is expressed as (32) Proof: Please refer to Appendix B.
Accordingly, problem (P1-dual) can be further modified as shown below.
Subject to (28), Remark 4.1: Note that some useful insights can be drawn from the results in Proposition 4.1.
First, with the dual variables given,f (D khk /λ 0 ) andf (A kḡk /λ k ) can be, respectively, interpreted (in terms of the dual function (27)) as the optimum offloading rate to helper k and the optimum results downloading rate from helper k, while 3 B k /(2λ k κ k ) and 3 (β 0 − ζ 0 )/(2λ 0 κ 0 ) represent the optimum computation frequencies at the kth helper and the local user, respectively.
Accordingly, when helper k enjoys good offloading (downloading) channel gain, the optimum offloading (downloading ) ratef (D khk /λ 0 ) (f (A kḡk /λ k )) also gets large due to non-decreasing monotonicity of W 0 (x). Moreover, provided that the total energy constraint for the local user is violated, thereby incurring a larger Lagrangian multiplier λ 0 (c.f. (21)), the optimum offloading time (t of f k 's) and the optimum local computation time (t c 0 ) under the same Π turn out to be longer. Hence, the total energy consumption for the local user gets reduced, which complies with Lemma A.1.
In a sum, given an initial (feasible) set of dual variables, the optimal solution to (27) is first obtained leveraging Proposition 4.1, and then the dual variables are readily updated utilizing some sub-gradient based method, e.g., ellipsoid method [28], until a predefined threshold controlling the accuracy of the algorithm is achieved.

B. Primal-Optimal Solution to (P1)
We aim for solving (P1) in this subsection. Since there might exist multiple solutions to (LP1) in each iteration of the ellipsoid method while there is only one optimal solution to the convex problem (P1), we retrieve the primal-optimal Π to (P1) through the dual-optimal solution in the sequel. Denoting the optimum dual variables by (η * , β * 0 , λ * 0 , ζ * 0 , λ * T , β * T , ζ * T ), the primal-optimal solution are related to the dual-optimal one as follows (c.f. (29) and (30)).
where A * k 's, B * k 's, and D * k 's are obtained by plugging the optimum dual variables into (25)  Note that (18c) and (18d) should be satisfied with 3 respectively, independent of Π, and thus can be safely removed from the constraints of (LP2).
Denoting the optimal solution to (LP2) by Π * , the primal-optimal solution to (P1) are thus given by

C. Sub-Optimal Solution to (P0-Eqv)
In this section, we propose a suboptimal scheme to jointly optimize task assignment as well as time and power allocation for (P0-Eqv) based on the optimal solution to (P1) developed in the previous subsections.
First, we propose to round off π * (l, k)'s as follows such that (18g) are satisfied: wherek l = arg max k∈K∪{K+1} π * (l, k) To ensure that each helper is assigned at least one task, we need to further adjust Π * to avoid the cases where some WD is assigned with no task after rounding off as (38). The detailed procedure of constructing such Π is shown in Fig. 2.
Next, given the updated Π * , what remains to be solved reduces to problem (P2). Hence, it can be solved using Lagrangian dual decomposition following similar procedures as shown in Section IV. A. The proposed joint task assignment and wireless resource allocation scheme is thus summarized in Algorithm 1.

D. Complexity
The complexity of Algorithm 1 mainly lies in solving (P1). Hence, we focus on discussing the complexity of solving (P1). Specifically, the complexity of solving (P1) includes solving (P1-dual ′ ) by ellipsoid method via primal-dual iterations, and solving (LP1) by simplex method in each iteration of the ellipsoid method. In accordance with the worst-case number of iterations for the ellipsoid method and the expected complexity of the simplex algorithm 5 , the complexity 5 Despite that the worst-case complexity of simplex algorithm is known to be exponential, it was shown in [30] that most LP can be approximated with their inputs perturbed and then solved by simplex algorithm in polynomial time.
of Algorithm 1 can be estimated by where H max h∈∂g(z), z∈E (0) h is a Lipschitz constant for (27) over the initial ellipsoid over E (0) , and ǫ is a parameter controlling the accuracy of the algorithm.

V. LOW-COMPLEXITY BENCHMARK SCHEMES
In this section, we develop two low-complexity benchmark schemes, which are provided in Sections V-A and V-B, respectively.

A. Fixed-Frequency Task Assignment and Wireless Resource Allocation
In this subsection, we consider a benchmark scheme that alleviates the WDs from adjusting their computation frequencies using DVFS by endowing them with the maximum computation capacities, i.e., f 0 = f max 0 , and f k = f max k , k ∈ K. Note that with the computation frequencies fixed, this scheme reduces to a special case of the joint task assignment and wireless resource allocation scheme with the constraints (18c) and (18d) being active. Since we have studied a similar fixed-frequency design in [1], in the sequel, we only focus on some results that are distinguished from Section IV due to space limitation.
First, with the computation frequencies fixed, the computation energy consumed by the local user and the helper turn out to be E c respectively. Moreover, t c 0 and t c k 's in problem (P0) are safely removed by replacing them with the left-hand side (LHS) of (18c) and (18d), respectively. In addition, since t c 1 cannot be adjusted once helper 1 is assigned with its tasks, I 1 (c.f. (14)) cannot be further simplified. To sum up, problem (P0-Eqv) reduces to be: Remark 5.1: Compared with problem (P0-Eqv), problem (P0 ′ -Eqv) is more sceptical to infeasibility 6 , since the minimum of E c 0 (E c k 's) cannot reach zero even with t c 0 (t c k 's) approaching to infinity (c.f. (3) ((9))). Once it is feasible, the solution to (P0 ′ -Eqv) can be found by a simplified version of Algorithm 1 that yields an approximate complexity of O(2(2K + 3) 2 log( √ γH/ǫ)P(L(K + 1))), which is lower than (39) due to decrease in the number of dual variables.

B. Greedy Task Assignment Based Wireless Resource Allocation
The joint task assignment and wireless resource allocation and its special case with fixed computation frequencies both admit complexity involving a polynomial factor P(L(K + 1))due to solvingπ(l, k)'s using simplex method. To further reduce computation complexity, we propose in this subsection greedy task assignment based wireless resource allocation. The main idea of this heuristic scheme is as follows. First, we obtain the objective value of problem (P2) as a cost function of the initial task assignment. Then we assign one task each time from those unallocated to a WD who yields the least amount of increase in the cost function. Next, since two different task selection criteria are applied, we compare the results of these two sub-schemes and choose the one that achieves the lower total latency.
The greedy task assignment based wireless resource allocation scheme is shown in Algorithm 2. In Algorithm 2, it is seen that for the first K tasks, each of them is in turn assigned to the helper with the best channel condition among those who have not yet been occupied. The intuition behind such assignment is that the selected helper will consume the least amount of energy in transmission (c.f. (7) and (13)), and thus any spare energy can be exploited for further latency reduction. It is also worth noting that the task with the longest (input/output) data flow is executed locally for the sake of saving data-transmission time.
The complexity of Algorithm 2 comprises that of two sub-schemes (with sorting {T l } and {R l }, respectively), each of which requires solving problem (P2) for (L−K −1)(K +1) times to

VI. NUMERICAL RESULTS
In this section, we provide numerical results to validate the effectiveness of the proposed joint task assignment and wireless resource allocation (Joint optimization) in Section IV, as compared against the fixed-frequency scheme (Fixed frequency) and the greedy task assignment based algorithm (Greedy assignment) presented in Section V as well as other benchmark schemes as follows.
• Optimal The MINLP problem (P0-Eqv) is solved by exhaustive search over feasible Π for 7 with a convex problem (P2) solved each time for a given Π. Note that "Optimal" is of exponential complexity and is thus too costly to implement in practice. Hence, we only provide this scheme for one numerical example in the sequel.
• Random assignment In this scheme, a random matrix with its entries drawn from i.i.d.
uniform distribution over [0, 1] is first generated, and then the procedure shown in Fig. 2 is employed to construct a feasible Π to (P0-Eqv). Next, solve problem (P2) under this Π.
• Local execution All of the computation tasks in L are executed locally with the total latency expressed as max by combining (18a) and (18c).
In simulations, the K helpers are located with their distance uniformly distributed over We consider a simple scenario where the local user has L = 5 tasks to be executed in the present of K = 2 helpers. Fig. 3 shows the total latency versus the energy constraints at the helpers assuming E 1 = E 2 . It is observed that the optimal scheme outperforms all the other ones, while our proposed "Joint optimization" achieves the second lowest total latency with little gap to the optimal solution. "Greedy assignment" admits reducing total latency with larger helpers' energy constraints, and outperforms "Local execution" in most cases except for those with E k 's below around −38dB, since when the helpers also suffer from a scarcity of energy, full local execution with constant total latency of max is intuitively better than computation offloading. Moreover, "Fixed frequency" does not work until E k 's is larger than −32dB due to the infeasibility of (P0-Eqv). Considering the same scenario as in Fig. 3, Fig. 4 shows the total latency versus the energy constraint at the local user with E 1 = E 2 = −10dB. It is seen that the average total latency decreases over E 0 for all of the schemes. In particular, "Joint optimization", "Greedy assignment" and "Fixed frequency" admit sharp reduce in the average total latency when E 0 is less than about −34dB, and slightly go down when E 0 continues decreasing. This complies with intuition as follows. When the energy resource is scarce at the local user, a little more energy will cause significant decrease in the computation offloading time, while as for E 0 beyond −34dB, the bottleneck mainly lies in the helper' energy constraint E k 's. It is also worth noting that "Fixed frequency" is substantially outperformed by "Joint optimization" and "Greedy assignment" when E 0 is less than about −34dB, since local computing with its full computation capacity is far from optimality under circumstances of limited energy supply. In addition, for "Local execution", l∈L C l f max 0 starts taking effect when E 0 is larger than about −31.5dB. outperform "Greedy assignment". This is because when f max k 's is below 1.32GHz, the helpers' optimal computation frequencies obtained by "Joint optimization" tend to be their maximum capacities, i.e., f max k 's, which is almost equivalent to "Fixed frequency" (c.f. Section V-A).
However, when f max k 's continuously increases, "Joint optimization" will suppress the helpers' computation capacities strictly below their limits such that under the given helper's energy budget, there is still sufficient energy left for results downloading thus achieving the overall better performance. In addition, "Local execution" under this setting is equal to 39.5ms, which is too large to illustrate and is thus removed from the figure herein.   The trade-offs between the total latency and the bit-length of the task input/output data are demonstrated in Fig. 6 with different energy constraints under the setup of equal input/output length across all tasks, i.e., T 1 = T and R 1 = R, ∀l ∈ L. It is observed that the schemes of "Joint optimization", "Fixed frequency", and "Greedy assignment" almost overlap each other. Among these three schemes, "Greedy assignment" is slightly worse than the other two in the lower range of data length in Fig. (6a), while "Fixed frequency" is noticeably worse than the other two in the higher range of data length in Fig. (6b). "Random assignment" continues being outperformed by the above three task offloading schemes especially when T (R) gets larger, which validates the importance of effective task assignment. It is also worth noting that with the energy setting of E 0 = −30dB at the local user and E k = −10dB, ∀k ∈ K, in the helpers, when T (R)

B. The Effect of Computation Load on the Average Total Latency
is larger than around 5.5 × 10 3 bits, "Local execution" becomes favourable due to increasing energy consumption in data transmission. By contrast, when the energy budget of the helpers increases to 0dB in Fig. (6b), all the task offloading schemes considerably outperform "Local execution", since with sufficient energy supply (c.f. (18b)), it costs the helpers less amount of time to transmit even very long task-output data.  The effect of the amount of computation required per task on the total latency is shown in Fig. 7 with different energy constraints' setup assuming C l = C, ∀l ∈ L. 8 It is seen that 8 Due to the relatively high regime into which the total latency of "Local execution" falls in Fig. 7, we present the results in tables to improve readability.
"Fixed frequency" turns out to be unavailable when C is larger than about 7.4 × 10 6 cycles and 6.0 × 10 6 cycles in Figs. (7a) and (7b), respectively, due to infeasibility incurred for the same reason as discussed in Fig. 3. It is also observed that the schemes of "Joint optimization", "Greedy assignment" and "Random assignment" increase with C very slowly when it is below about 5.0 × 10 6 cycles (4.0 × 10 6 cycles) in Fig. (7a) (Fig. (7b)), which is because under low to medium computation load per task, very little increase in t c 1 is already sufficient to satisfy (19b). However, when C gets larger, the increasing amount of computation energy yields decrease in communications energy thus prolonging the total latency. In addition, "Local execution" admits the best performance among all the schemes when the computation burden is below 3.1 × 10 6 cycles per task under the lower energy setting of the helpers in Fig. (7a), since "Local execution" is preferred for computation non-intensive tasks. However, when there is drastic difference in the energy budget between the local user and the helpers, the advantage of the task assignment schemes is prominently observed in Fig. (7b), as similarly explained for Fig. (6b).  that while the average total latency achieved by "Joint optimization", "Greedy assignment", and "Fixed frequency" steadily increase with the total number of tasks, that of "Local execution" grows fast, which motivates cooperative MEC especially when there are a large number of tasks to be executed. Moreover, the performance of "Random assignment", almost worst among all the schemes, stresses the importance of proper task distribution schemes. "Greedy assignment" is seen to strike satisfying balance between the performance and the complexity.

VII. CONCLUDING REMARKS
In this paper, we investigated the joint task assignment, communications rate, as well as computation frequency allocation for a D2D-enabled multi-helper MEC system assuming binary task offloading. Under a TDMA communication protocol, we aimed for minimizing the overall computation latency subject to individual energy and computation capacity constraints at both the local user and the helpers. Since the formulated problem was an MINLP that is in general difficult to solve, we proposed an efficient convex-relaxation-based algorithm to construct a suboptimal task assignment solution based on the optimal one to the relaxed problem. Furthermore, a benchmark scheme with fixed computation frequency and a greedy task assignment based heuristic algorithm were also developed to strike the balance between complexity and performance. Finally, numerical results verified that our proposed design is an effective solution to enhance the local user's computation latency, by exploiting D2D collaborations at the network edge.
Due to space limitation, there are several other challenging issues yet to be addressed in this paper, which will be investigated in our future work. First, this paper assumed that all of the WDs are deployed at fixed locations with static wireless channels. In practice, since these WDs may move over time, the wireless channels will fluctuate over time, and the D2D connections established between the local user and the helpers may also be dropped as they move away from each other. Under these circumstances, we need to consider new design principles (e.g., online algorithms with long-term objectives capturing dynamics of the D2D links) to combat such mobility issues [24]. Second, in this paper we considered TDMA owing to its easy implementation in practice. Other orthogonal multiple access schemes, e.g., OFDMA [13,22], and more sophisticated NOMA schemes, e.g., sparse code multiple access (SCMA) [31], can be employed to further enhance the system performance. Under these schemes, on top of mixedinteger task assignment, how to design subcarrier allocation for OFDMA and joint message decoding for NOMA are also quite challenging issues worthy of further study. In addition, our considered model assumed that the helpers have agreed to cooperate in the computation offloading, while we believe there are also different types of incentive driven collaboration schemes that require sophisticated design. At last, it is also worth investigating how to extend our current single-user multi-helper MEC model to a multi-user multi-helper one by proper design of multi-user scheduling algorithms. First, in order to prove Proposition 3.1, we need the following lemma.
Lemma A.1: The function h(y, t) = f y t t monotonically decreases over t > 0. Proof: The monotonicity can be obtained by evaluating the first-order partial derivative of h(y, t) w.r.t. t, and using the fact that (1 − x)e x − 1 < 0, for x > 0.
On one hand, there are two possible cases for the optimal I k 's given by (15): case 1) k j=1 t of f j + t c k > I k−1 + t dl k−1 ; and case 2) k j=1 t of f j + t c k ≤ I k−1 + t dl k−1 . In line with Lemma A.1, the total transmitting energy of the kth helper, i.e., E dl k 's (c.f. (13)), ∀k ∈ K, monotonically decreases over t dl k 's. Hence, if the first case occurs, helper k − 1 (k ≥ 2) can slow down its downloading, i.e., extending t dl k−1 , until k j=1 t of f j +t c k = I k−1 +t dl k−1 is satisfied (c.f. Fig. 1), such that I k remains unchanged but the transmitting energy of helper k − 1 gets reduced. As such, the two cases can be, w.l.o.g., merged into one as k j=1 t of f j + t c k ≤ I k−1 + t dl k−1 , ∀k ∈ K \ {1}, which suggests the following constraint on t c k given by As it is seen clear that (15) reduces to the waiting time for helper k (k ≥ 2) can be recursively obtained as On the other hand, there are also two possible cases for the optimal I 1 given by (14): case 1) t of f 1 + t c 1 < K k=1 t of f k ; and case 2) t of f 1 + t c 1 ≥ K k=1 t of f k . Since the kth helper's computation energy E c k given by (9), ∀k ∈ K, monotonically decreases over t c k 's, if case 1) takes place, it is always possible for helper 1 to slow down its computation such that t of f 1 + t c 1 = K k=1 t of f k is met without violating the other constraints. In a sum, I 1 , w.l.o.g., reduces to subject to (19a). Combining (44) with (45), and substituting the results for I k−1 in (42), (19b) follows.
Furthermore, substituting (45) for I 1 in (44) with k = K, it is easy to obtain that I K = t of f 1 + t c 1 + K−1 j=1 t dl j , and thus the completion time T reduces to T = t of f 1 + t c 1 + j t dl j . As a result, the total latency given by (17) turns out to be T total = max{t c 0 , t of f In addition, it can also be verified that when the optimal T total given by (46) yields t c 0 > t of f 1 + t c 1 + K k=1 t dl k , it is always possible for one of the K helpers to slow down its downloading rate such that t c 0 = t of f 1 + t c 1 + K k=1 t dl k without violating the other constraints. Therefore, T total , w.l.o.g., can be further simplified as Since the solution to the equation f (x) − xf ′ (x) = y for x > 0 is shown to be x =f (−y) [27], it is easy to verify that (48a) and (48b) imply (29) withπ(l, k) replaced byπ(l, k) therein.
Next, to find an optimum solution to (27) in terms of Π we substitute for t of f k , t dl k , and t c k , respectively, in (24), such that (22) is expressed in terms ofπ(l, k)'s. Furthermore, consideringπ(l, k)'s as the only (primal) variables, after some manipulations, the minimization of (22) subject to (18e), (18f), and (20) is formulated as problem (LP1), which can be solved using simplex method. Denoting the optimal solution to (LP1) asπ(l, k)'s, the results given by (29) and (30) are thus obtained, which completes the proof.