extfieldsofstudy:
Mathematics
Computer Science

s2fieldsofstudy:
Computer Science
Mathematics

text:
Domain Decomposition Methods for Space Fractional Partial Differential Equations

In this paper, a two-level additive Schwarz preconditioner is proposed for solving the algebraic systems resulting from the finite element approximations of space fractional partial differential equations (SFPDEs). It is shown that the condition number of the preconditioned system is bounded by C(1+H/\delta), where H is the maximum diameter of subdomains and \delta is the overlap size among the subdomains. Numerical results are given to support our theoretical findings.

with p i , q i ∈ R satisfying p i = q i and p i + q i = 1. We may find that equation (1.3) can be obtained from (1.1) by takingM (1.4) where e i is the ith column of identity matrix in R d×d and δ denotes the Dirac function on S d−1 . Extensive numerical methods have already been developed for SFPDEs, like finite difference methods (see e.g., [2,17,18,25,26]), finite element methods (see e.g., [6-8, 14, 22]), and spectral methods [11]. Due to the nonlocal properties of fractional differential operators, the most important issue for numerical computation of SFPDEs is how to reduce the computation costs. Some methods for the reduction have already been designed, like alternating-direction implicit methods (ADI) [18,30,31], special iterative methods [12,20,[31][32][33][34] and multigrid methods [20,36].
The discrete systems Ax = b of SFPDEs usually have the following characteristics: 1). the condition number of A increases fast, as the mesh becomes fine; 2). the coefficient matrix A is dense. As for iterative methods, two issues need to be concerned for efficiency: one is how to construct good preconditioners for the discrete system Ax = b, which may help us to save the iterative steps; the other is how to reduce the computation cost of each iterative step, for which some papers employ the multiplication of Toeplitz matrices and vectors with n log(n) computation complexity (see e.g., [12,20,29,[31][32][33][34]). If the iterative step can be carried out in parallel, the efficiency of solving SFPDEs may be significantly improved. The parallelizable algorithms have been wildly used in numerical solutions for PDEs (see e.g. [28]). Whereas, to the best of our knowledge, no parallelizable algorithms have been designed for SFPDEs.
In this paper, we shall construct a two-level additive preconditioner for the discrete system Ax = b resulting from the finite element approximation of (1.1), and then use the preconditioned conjugate method (PCG) to solve it. The preconditioner we construct is almost optimal, i.e., the condition number of the preconditioned system is bounded by C(1 + H/δ), where H is the maximum diameter of subdomains and δ is the overlap size among the subdomains. Moreover, the preconditioner may be employed in parallel or each step of the PCG can be carried out in parallel. As a result, the whole numerical solution processes, including the generation of A, b and the multiplication of matrices and vectors, may be conducted in parallel.
Without loss of generality, we focus on the case d = 2, namely, we consider the problem (1.1) in R 2 . For Λ ⊂ R 2 , denote L 2 (Λ) the space of all measurable function v on Λ satisfying and they are abbreviated as (v, w) and ||v|| respectively if Λ = R 2 .
To simplify our statement, we make a convention here: function v defined on a domain Λ ⊂ R 2 also denotes its extension on R 2 which extends v by zero outside Λ. The constant C with or without subscript shall denote a generic positive constant which may take on different values in different places. These constants shall always be independent of mesh sizes and numbers of subdomains. Following [35], we also use symbols , and ≈ in this paper. That a 1 b 1 , a 2 b 2 and a 3 ≈ b 3 mean that a 1 ≤ C 1 b 1 , a 2 ≥ C 2 b 2 and C 3 b 3 ≤ a 3 ≤ C ′ 3 b 3 for some positive constants C 1 , C 2 , C 3 and C ′ 3 . The rest of the paper is organized as follows: in section 2, the variational problem of (1.1) and its finite element discretization are described; in section 3, the two-level additive preconditioner for the SFPDEs is presented; in section 4, it is proved that the preconditioner is almost optimal; finally, in section 5, the numerical results shall be given to support our theoretical findings.

The model problem and its discretization
In this section, we shall describe the SFPDEs in details, and then introduce its variational formulation and finite element discretization.

The model problem
The µth order fractional integral in the direction z = (cos θ, sin θ) is defined by where Γ is the Gamma function.
[8] Let n be a positive integer, and θ ∈ R. The nth order derivative in the direction of z = (cos θ, sin θ) is given by Let n be the integer such that n − 1 ≤ µ < n, and define σ = n − µ. Then the µth order directional derivative in the direction of z = (cos θ, sin θ) is defined by If v is viewed as a function in x, D µ 0 , D µ π are just the left and the right Riemamm-Liouville derivatives (see Appendix).
For u : R 2 → R, define differential operator L α in R 2 as Denote Ω a polygonal domain in R 2 , set 1/2 < α < 1. The model problem of this paper is to find u :Ω → R such that L α u = f, in Ω, u = 0, on ∂Ω, where f is a source term and we assume thatM (θ) satisfiesM (θ) =M (θ + π) for θ ∈ R, i.e., (2.1) is a symmetric problem. Here, we recall the convection made in Section 1: u also denotes its extension by zero outside Ω.
The weak derivative D µ θ v is unique if it exists and the weak derivative coincides with the correspondent derivative defined in Definition 2.3 if v ∈ C ∞ 0 (R 2 ). In the following, we use D µ θ v to denote the weak derivative.
We construct a quasi-uniform triangulation Divide each Ω i into smaller simplices τ j of diameter O(h), such that Γ h = {τ j } form a finer triangulation of Ω. Denote V H and V h piecewise linear finite element function spaces defined on the triangulations Γ H and Γ h respectively. It is known that V h ⊂ H 1 0 (Ω) ⊂ H α 0 (Ω). The finite element approximation for (2.4) (the details please see [9]) is to find is a symmetric bilinear form, and But the discretization is hardly computed. So we use (2.5) instead of (2.7), where B(·, ·) is understood as the approximation toB(·, ·). For the details, please refer to [9].

A two-level additive Schwarz preconditioner
Since B(v, w) is a symmetric bilinear form, by (2.6), we know that A : V h → V h is symmetric positive definite with respect to (·, ·), i.e., Then bilinear form also induces an inner product on V h . Set norm The problem (2.5) can be restated as to find u h ∈ V h such that For the above equation, we shall construct our two-level Schwarz preconditioner and then use PCG method to solve it. Our preconditioner is designed by making use of the following overlapping domain decompo- where the subdomain Ω ′ i contains coarse subdomain Ω i , and satisfies the diameter of Ω ′ i ≈ O(H). Meanwhile the boundary of Ω ′ i align with the mesh of triangulation Γ h , and the distance from ∂Ω ′ i ∩ Ω to Ω i is greater than δ , which is a positive constant measuring the overlapping size among the subdomians. Define subspaces V i of V h as For our analysis, we regroup the subregions in terms of the coloring strategy (see e.g. [28]). By a minimal or good coloring, we group the subregions {Ω ′ i } into J C classes, each of which has some disjoint subregions and can be regarded as one subregion. Exactly, decompose the index set

It is not hard to verify that
To help our analysis, we define some projectorsP i : Remark 3.1. Different from the integer order PDEs, it is interesting to see thatP i = j∈I i P j . Now, we are ready to present our two-level additive Schwarz pre-conditioner, i.e.,

By (3.4), we have
Define P h := J i=0 P i , and then the preconditioned system is In the next section, we shall prove the condition number of P h is bounded by C(1 + H/δ), where the constant C is independent of mesh size and the numbers of subdomains, but dependent of J C .

Condition number estimate
We first introduce two interpolation norms and relevant Sobolev spaces (see e.g., [27]). Let Λ be a domain in R 2 . For integer m, denote by || · ||H m (Λ) the Sobolev norm of integer order m, i.e., Let µ > 0 be a non-integer and 0 < s < 1, n is a non-negative integer such that n < µ < n + 1. We introduce the interpolation norms

Relevant Sobolev spaces arẽ
For µ > 0, it is known thatH µ (R 2 ) coincides with H µ (R 2 ). The following norms relation is useful in our analysis. In the following, we shall give some useful results. Proof. Without loss of generality, we prove that (4.3) holds under the situation as the figure 1 shows. Cartesian coordinate x ′ Oy ′ is obtained by rotating xOy θ angle counterclockwise. We set that p i ∈Ω i and p j ∈Ω j such that the length of line segment p i p j is equal to dist(Ω ′ i , Ω ′ j ), the graphs of functions x ′ = Γ i1 (y ′ ) and x ′ = Γ i2 (y ′ ) are parts of boundary ∂Ω ′ i such that Ω ′ i = {(x ′ , y ′ ) : Γ i1 (y ′ ) < x ′ < Γ i2 (y ′ )}, the graphs of functions Γ j1 (y ′ ) and Γ j2 (y ′ ) are parts of boundary ∞ are the left and right Riemann-Liouville fractional derivative operator (see Appendix). We in the following only prove (4.3) for the case that there exists a ray in (cos θ, sin θ) direction going through both Ω ′ i and Ω ′ j from Ω ′ i to Ω ′ j . Indeed, when the other case is true, (D α θ v, D α θ+π w) = 0 and (4.3) naturally holds. Denote The shaded area is where in the fifth equality we have used the relation Γ(1 − α) = −αΓ(−α), the first inequality is by (x ′ − s) ≥ lH/2 when s ≤ Γ i2 (y ′ ) and the last by the Cauchy-Schwarz inequality. Then where the second inequality is by Denote d i as the diameter of Ω ′ i , and define a function in x ′ as where v * w denote the convolution product (see e.g., [1] where the second inequality is by Lemma 2.7, the last equality is by (3.2). Combining (4.9) with (4.6), we obtain (4.10) Similarly we may also obtain (4.11) By Lemma 2.7, (4.13) Proof. Note that it is easy to see that the lemma follows after we prove that holds for any θ. So next, we give a proof of (4.14). It is easy to see that In fact, we may prove through using the following two kind of inequalities: (by Lemma 2.7 and (3.2)); Fixing i ∈ S, we write It is easy to see that the sum of the terms containing v i on the left hand of (4.16) We know that card(S i ) ≤ C, (4.21) where card(S) denotes the number of elements contained in S, C is a positive constant dependent of J C . By (4.21) and (4.17), we have As figure 2 shows, L i , containing Ω ′ i , denotes the shaded bar area whose edges are parallel to the vector (cos(θ), sin(θ)) and go throughΩ ′ i . It is not hard to see that (D α H} is bounded by a positive C which is a positive constant dependent of J C . And then using (4.18), we have From above analysis, using (4.17) and (4.18) to enlarge the correspondent terms of the left hand of (4.16), we may show that the inequality (4.16) holds with constant C being independent of θ. ✷ Lemma 4.4. There exist linear operators κ i (i = 0, 1, . . . , J), κ 0 mapsH 1 0 (Ω) into V 0 , κ i (i = 1, . . . , J) mapsH 1 0 (Ω) and V h intoH 1 0 (Ω ′ i ) and V i respectively, such that for any v ∈H 1 where the constant C depends on J C .
Lemma 4.5. The operator P h is symmetric positive definite with respect to the inner product (·, ·) A .
Proof. For the proof, we refer to section 2.3 in [28].
Proof. By the definition of P 0 , we have We recall the index set I i (i = 1, 2, . . . , J C ): for any j, k ∈ I i with j = k, Ω ′ j ∩ Ω ′ k = ∅, so it is easy to see that For any w ∈Ṽ i , decompose it as w = j∈I i w j , w j ∈ V j . Then by (4.37), we have Taking w j = P j v in the above equation, and using Cauchy-Schwarz inequality and Lemma 4.3, we may obtain It is easy to see that Combining (4.42) with (4.36), we arrive at Then the lemma is proved. ✷ Proof. Let κ i be the operators as defined in Lemma 4.4, we now prove the decomposition Actually, using the definition of the norms || · ||H α (R 2 ) , || · ||Ĥ α (Ω) , Lemma 4.1 and Lemma 4.4, we have Proof. By Lemma 4.7, we know that there exists a decomposition v = J i=0 v i , v i ∈ V i such that (4.44) holds. Then by the Cauchy-Schwarz inequality and (4.44), we have which is Lemma 4.8 ✷.
Theorem 4.9. The condition number of P h satisfies Proof. The theorem follows from Lemma 4.6 and Lemma 4.8. ✷

Numerical results
We test PCG method with B h as a preconditioner. The tests are carried out by using Matlab software. The stopping criterion is where u k is obtained from u k−1 by one step PCG iteration, and u 0 = 0. Here we only test the case that Ω is a square domain. We take the triangulations Γ H , Γ h as uniform ones as in figure 3 shows. With viewing the one in figure 3 as the coarse triangulation Γ H , the subdomains Ω i and Ω ′ i are taken as the darker and lighter shaded squares respectively. We give two examples in this paper: one is with the probability measureM having a discrete form and the other is withM being a continuous function. Table 1 and Table 2 list the iterative steps of our PCG methods for Example 5.1 and Example 5.2 respectively. It can be seen that when we fix H h = constant, for instance H h = 8, the convergence rate of our PCG method is optimal, i.e., the iterative steps of the PCG method almost keep to be a constant, which coincides with our theoretical results in this paper. From the tables 1-2, we also find that small overlapping case is a good choice for our domain decomposition methods for the SFPDEs, which is same as the overlapping domain decomposition methods for the integer order PDEs.
where f (x, y) is taken such that the solution is u = ((2x − x 2 )(2y − y 2 )) 4 , where n is an integer such that n − 1 ≤ µ < n. For p > 0, k is an integer satisfying k − 1 ≤ p < k.


extfieldsofstudy:
Computer Science
Mathematics

text:
Energy-Efficient Coordinated Multi-Cell Multigroup Multicast Beamforming with Antenna Selection

This paper studies energy-efficient coordinated beamforming in multi-cell multi-user multigroup multicast multiple-input single-output systems. We aim at maximizing the network energy efficiency by taking into account the fact that some of the radio frequency chains can be switched off in order to save power. We consider the antenna specific maximum power constraints to avoid non-linear distortion in power amplifiers and user-specific quality of service (QoS) constraints to guarantee a certain QoS levels. We first introduce binary antenna selection variables and use the perspective formulation to model the relation between them and the beamformers. Subsequently, we propose a new formulation which reduces the feasible set of the continuous relaxation, resulting in better performance compared to the original perspective formulation based problem. However, the resulting optimization problem is a mixed-Boolean non-convex fractional program, which is difficult to solve. We follow the standard continuous relaxation of the binary antenna selection variables, and then reformulate the problem such that it is amendable to successive convex approximation. Thereby, solving the continuous relaxation mostly results in near-binary solution. To recover the binary variables from the continuous relaxation, we switch off all the antennas for which the continuous values are smaller than a small threshold. Numerical results illustrate the superior convergence result and significant achievable gains in terms of energy efficiency with the proposed algorithm.

I. INTRODUCTION
The energy efficiency (EE) of the multi-antenna technologies is a performance indicator growing its importance for wireless communications and future 5G network design [1]- [6]. The total power consumption consists of various elements, such as base stations (BSs) direct current (DC) power, baseband processing, radio frequency (RF) processing, coding, decoding, and backhaul operations. In fact, the data transmit power takes only a small share of the total power consumption when moving towards 5G networks [7]. Each active antenna element requires an RF chain to preprocess the signal which has a significant impact on the power consumption. An efficient method to improve the energy efficiency is to use only some of the antennas for transmission by switching off a portion of the RF chains, giving rise to the problem of joint beamforming and antenna selection [5], [8].
Coordinated beamforming [9] has been one of the most efficient methods to deal with the inter-cell interference in the research on the multi-cell networks. The technique, adopted in the current LTE-A standards, uses the coordination of channel state information or inter-cell interference conditions to improve the system performance. Energy-efficient coordinated beamforming has been extensively studied recently, e.g., in [6], [10]- [12].
The research presented before focuses mostly on unicast beamforming where each user is assigned with an independent data stream. However, the increasing demand for high data rate applications such as video broadcasting services creates new challenges. This gives rise to the concept of multicast beamforming, where multiple users desire to receive the same information. Multicasting is a particularly powerful technique in the context of cache-enabled cloud radio access networks proposed for 5G systems, where it can be used to transmit the same popular contents to multiple users to improve both spectral and energy efficiency [13]. The physical layer multicasting has been also included in the LTE standards and it has applications in satellite communications [14], [15] as well. Multicast beamforming problems have been studied in single-cell systems for different optimization targets, e.g., transmit power minimization [16], [17], max-min fairness [17], [18], and sum rate maximization [19]. Joint beamforming and antenna selection for transmit power minimization was studied in [20]. Coordinated multicast beamforming for transmit power minimization and max-min fairness has been studied in [21]. He et al. [22] proposed a centralized algorithm for the energy efficiency maximization problem in multi-cell system with single-group per cell. However, [22] only considered the beamforming problem without taking into account the fact that significant energy efficiency savings can be achieved by switching off some of the RF chains, i.e., antenna selection.
In this paper, we study energy-efficient coordinated beamforming in multi-cell multigroup multi-user multicast multipleinput single-output (MISO) systems with joint beamforming and antenna selection. We first introduce a novel way to reduce the feasible set of the continuous relaxation, resulting in better performance compared to the original formulation. Then, by proper equivalent transformations, the continuous relaxation of binary variables, and the Charnes-Cooper transformation, we propose a successive convex approximation (SCA) based algorithm to solve the problem. Thanks to the proposed formulation, the continuous relaxation yields some continuous values being close to binary ones. To this end, we switch off antennas which have continuous values close to zero.
Contrary to [22], where an energy-efficient coordinated beamforming method for fixed antenna sets was derived with single group per cell, we focus on energy-efficient joint beamforming and antenna selection problem with multiple groups per cell which has not been considered in the previous multicasting research. In addition, we set antenna specific maximum power constraints whereas [22] only considered the BS-specific sum power constraints which can result in large dynamic range of the PA output powers, causing large phase deviations due to phase non-linearities [23].
The rest of the paper is organized as follows. Section II presents the system model, power consumption model and the considered optimization problem. The proposed algorithms are provided in Section III while numerical evaluation and conclusions are presented in Section IV and V, respectively.
The following notations are used in this paper. |x| denotes the cardinality of x if x is a set, and absolute value of x, otherwise. ||x|| 2 is the Euclidean norm of x and boldcase letters are vectors. x T , x H , Re(x) and Im(x) mean transpose, Hermitian transpose, real part and imaginary part of x, respectively.

A. System Model
We consider a multi-cell multigroup multicasting system with B BSs, where each BS b ∈ B = {1, . . . , B}, transmits independent information to G b groups of users in its cell. The set of groups served by BS b is denoted by G b = {j, j+1, . . . , j+G b −1}. The total number of users and groups in the network is denoted by K = |K|(K = {1, . . . , K}) and G = |G|(G = {1, . . . , G}), respectively. BS b is equipped with N b transmit antennas, whereas each user has only one receive antenna. The serving BS of user group g is denoted as b g . The set of users in group g is denoted by K g ⊂ K.
Since each user belongs to only one group, the sets of users belonging to different groups are disjoint, i.e., K i ∩ K j = ∅, ∀i, j ∈ G, i = j. The received signal at user k in group g is given by where h b,k ∈ C 1×N b is the channel (row) vector from BS b to user k, w g ∈ C N b ×1 is the transmit beamforming vector of group g, s g ∈ C is the corresponding normalized data symbol and n k ∼ CN (0, σ 2 k ) is the complex white Gaussian noise sample with zero mean and variance σ 2 k . The SINR of user k can be written as and the data rate towards user k is given as

B. Power Consumption Model
The total power consumption at the transmitter side is [5] where a b,i ∈ {0, 1} is the antenna selection variable for the ith transmit antenna of BS b, i.e., a b,i = 1 if the ith antenna is selected and a b,i = 0 otherwise, and η ∈ [0, 1] is the power amplifier efficiency. Note that the first term in (4) is the power consumption of the PA's to get the desired output powers, P RF is the power consumption of all circuit blocks in each active RF chain, P sta is the static power spent by cooling systems, power supplies, etc.

C. Problem Formulation
The optimization target is to maximize network energy efficiency whereΓ k is the SINR threshold for user k, includes the beamforming coefficients related to antenna i of BS b and index j refers to a group index served by BS b. In the numerator of the objective (5a), the rate for user group g is defined as a minimum of the user rates in the group, because common information is transmitted to the group. The constraint in (5b) guarantees that if a b,i = 0, then the beamformers associated with antenna i are set to zero. It also sets the maximum antenna specific power constraint. The above problem is a non-convex mixed-Boolean fractional program which is hard to tackle as such. One of the main challenges is that the problem is non-convex even when the binary variables are relaxed to continuous. More specifically, in that case, (5c) and the numerator of the objective function are non-convex.

A. Equivalent Transformation
To find a more tractable reformulation, we first equivalently transform (5) as max w,γ,v,a,r {r g } g∈G are new variables representing the SINR of each user k, the power transmitted from antenna i of BS b [5], and the rate of the weakest link in user group g, respectively. It is not difficult to see that (5) and (6) are equivalent, because the constraints in (6e) and (6b) are active and r g = min k∈Kg log(1 + γ k ), ∀g ∈ G in (6g) is satisfied at the optimum. We have equivalently replaced r g ≤ min k∈Kg log(1+Γ k (w)) with r g ≤ log(1+Γ k (w)), ∀k ∈ K g in (6g), because if the minimum user rate satisfies the constraint, then all the user rates in that group have to satisfy it. In constraint (6b), we have introduced a new formulation can be interpreted as a penalty parameter which penalizes the values of a b,i so that they are encouraged towards binary solution in case of continuous relaxation. The new formulation (6b) reduces the feasible set of the problem when continuous relaxation is performed, so that the larger α, the smaller feasible set. In other words, we have the following inequalities where EE bin , EE cont,α > 1 , EE cont,α = 1 and EE cont,orig are the optimal objective of Boolean formulation (5), continuous relaxation of (6) with α > 1, continuous relaxation of (6) with α = 1 and continuous relaxation of (5). In fact, the proposed formulation is a generalization of the approach used in [5] for single-cell multiuser MISO unicasting, where α = 1 and α = 2 were used. However, the numerical results show that the value of α has a significant effect on the algorithm performance, and, thus, the energy efficiency.

B. Solving the Equivalent Transformation
By taking a look at (6), we observe that if the binary variables are relaxed as 0 ≤ a b,i ≤ 1, the objective function is a concave-convex fractional function and the difficulty in solving (6) is due to the constraints (6e) and (6b). By using the continuous relaxation and introducing new variables β {β k } k∈K to represent the total interference-plus-noise of user k, the problem becomes max w,γ,v,a,β,r Now the left side of (8d) is linear and the right side is convex quadratic-over-linear function. To deal with the right side of (8d), we can write its linear lower bound approximation at point (w To find a more tractable formulation for (8b), we first equivalently write it as where the left side is a convex quadratic-over-linear function and the right side is convex. We can write the linear lower bound approximation of the right side at point a With the approximations (9) and (11), the problem becomes max w,γ,v,a,β,r At this point, we note that the problem is a concave-convex fractional program which can be transformed to a convex one with the Charnes-Cooper transformation [24]. Thus, solving (12) is equivalent to the following convex problem max φ,w,γ,v,ā,β,r g∈Gr g (13a) After solving (13), the optimal solutions for the original fractional program (12) can be found as w * g are the optimal variables of (13). As a result, we use the idea of successive convex approximation [25], where the convex problem (13) is solved at iteration n. After solving the problem at iteration n, the optimal solutions w * g , β * k , a b,i are then used to update Ψ (n+1) k (w g , β k ) and χ (n+1) b,i (a b,i ) for the next iteration. The monotonic convergence of the objective function (13a) is not difficult to see, and a more detailed convergence analysis for the problem with similar structure can be found, e.g., in [26, Appendix A]. Remark 1. If at least one of the SINR targetsΓ k in some user group of BS b is non-zero, we can further reduce the feasible set of (13) by adding the constraints where X b is the number of groups served by BS b which have at least one user having non-zero SINR target. This can be done because it is known that at least X b antennas have to be active to be able to transmit X b independent data streams.
Recovering the Binary Solution from Continuous Relaxation: Generally, solving the continuous relaxation usually results in a solution where many of the antenna selection variables are non-binary. However, due to the new formulation in (8b), many of the continuous antenna selection variables converge to binary. Thus, we propose to switch off all the antennas for which a b,i < ǫ, where ǫ is a small threshold. After performing the antenna selection, the algorithm needs to be run again with the selected antenna set to find the beamformers with lower dimensions. However, in the numerical results, we also illustrate that the beamformers returned by the relaxed problem already yields a good energy efficiency with the correct choice of α. This method is called 'Alg. 1 'simple'' in the numerical results. The proposed joint beamforming and antenna selection method is summarized in Algorithm 1.
Initial Points: One of the challenges is to find feasible initial points especially due to the QoS constraints. The initial a (0) can be set to all-ones. To find w (0) , β (0) , we can apply a similar approach as that in [27]. We first initialize any w (0) , β (0) and then consider a penalized formulation of (13), which can be written in a general form as where f 0 (x) is a concave objective function, f i (x) and F j (x, x (n) ) are convex constraints, x includes all the variables in (13), λ is a positive penalty parameter and q ∈ R Q1+Q2 + are new slack variables, one for each constraint in (13). The above problem is iteratively solved by updating the constraints F j (x, x (n) ) (i.e., Ψ (n) (13)) after each iteration, until q i = 0, ∀i. The feasible point is found very efficiently Algorithm 1 Proposed joint beamforming and antenna selection design.

4:
n := n + 1. 5: until convergence Output: a * b,i =ā * b,i φ * , ∀b ∈ B, i ∈ N b 6: Set a b,i = 0, for all b, i for which a * b,i < ǫ. 7: Run steps 1-5 again with fixed a to find beamformers with reduced dimensions. Output: w * g =w * g φ * , ∀g ∈ G because the elements of q are encouraged to zero due to the penalty function in the objective.

IV. NUMERICAL RESULTS
We evaluate the performance for a quasistatic frequency flat Rayleigh fading channel model with B = 2 base stations. Each BS serves two groups of users with |K g | = L = 2 users per group, i.e., the total number of users in the network is 8. The worst-case interference scenario is assumed so that the average path loss from all the BSs to all the users is 0 dB. We assume a unit bandwidth and noise power is normalized to σ 2 k = 0 dBW, ∀k ∈ K. We set N b = N for all b, i.e., N is the number of antennas at each BS, and the algorithms have been stopped when the change of the objective value has been smaller than ξ = 10 −6 between two last iterations. We set η = 0.35, ǫ = 0.001, P max = 9 dBW, P sta = 2 W, Γ k =Γ, ∀k ∈ K, and the other simulation parameters are given in the figures. Fig. 1 illustrates the convergence of the relaxed problem and the achieved energy efficiency for two different values of α, and two random channel realizations (denoted as 'CH' in the figure). The same initial points have been used for both values of α. We can see that the convergence speed is fast in the considered settings for both cases. However, we observe that they converge to different solutions in both examples. Specifically, the objective value of the relaxed problem after convergence is higher for α = 1, but the achieved energy efficiency (after recovering the binary solution) is worse than with α = 1.3. This is because with α = 1, more antenna selection variables are non-binary after convergence, which results in worse antenna selection result. Another observation is that with α = 1.3 (which achieves better energy efficiency), the objective value returned by the relaxed problem and the achieved energy efficiency are very close to each other. This means that the solution of the relaxed problem is already very close to binary. The examples demonstrate the effectiveness of the proposed formulation. The impact of α on the average energy efficiency is studied in the next experiment.   Fig. 2 illustrates the effect of α on the average energy efficiency forΓ = 0 dB andΓ = 6 dB. We also illustrate the performance of the simplified algorithm (Alg. 1 'simple'), where the beamformers achieved from the relaxed problem (i.e., step 7 is ignored in Alg. 1) are used for transmission. We can see that the choice of α affects the achieved energy efficiency, and the choice of the best α depends on the system parameters. In this case, α = 1.6 gives the best performance whenΓ = 0 dB, and for largerΓ, it is better to use larger α. The reason for this is that with largerΓ, the transmit (TX) power has to be higher to satisfy the SINR constraints. As a result, in the relaxed problem, also the relaxed variables a b,i have to be larger, implying that the larger 'penalty' (i.e. α) is required to force them to zero. We can also see that the larger the α, the better results are achieved with the simple method, implying that the relaxed problem yields near-binary solutions. This again demonstrates the benefit of the proposed formulation. Fig. 3 illustrates the average energy efficiency versus P RF . Specifically, the proposed algorithm is run with α = 1 and α = 1.5, and compared with the method where only beamformers are optimized (Alg. 1, no AS). First, we can observe that Alg. 1 with α = 1.5 achieves the best results, and the gap between α = 1 and α = 1.5 increases with P RF . It is interesting to observe that with α = 1.5, even  the beamformers obtained from the relaxed problem result in very good performance, which again verifies that the relaxed problem yields near-binary solutions. However, the simple method with α = 1 is even worse than the method without antenna selection for larger P RF . Finally, we can see that the proposed algorithm provides significant energy efficiency gains over the method without antenna selection (i.e., roughly 5-50% with the considered system parameters). Fig. 4 plots the average EE versus the number of antennas per BS. First, we can see that EE starts to decrease without antenna selection when N > 14. As can be observed, significant gains are achieved with the proposed algorithm and the gains increase with the number of antennas. Once again, significant gains are achieved with α = 1.5 compared to the method with α = 1. In addition, the simple method with α = 1.5 gives even better performance than the original algorithm with α = 1 and the gain increases with N . The simple method with α = 1, on the other hand, yields poor performance. Fig. 5 illustrates the transmit powers versus N with the proposed algorithm and a method without antenna selection. We can see that it is energy-efficient to increase sum power when the number of antennas increases. However, the average transmit power per active antenna decreases with increasing N also with antenna selection. This means that the more available antennas, the more antennas are chosen to maximize the energy efficiency.

V. CONCLUSIONS
This paper has studied energy-efficient multi-cell multigroup coordinated joint beamforming and antenna selection with antenna-specific maximum power constraints and userspecific QoS constraints. The resulting mixed-Boolean nonconvex optimization problem was tackled by continuous relaxation and successive convex approximation, where the antennas for which continuous antenna selection variables converge to zero are switched off. The numerical results have illustrated the achieved energy efficiency gains of the proposed methods over the method without antenna selection.


DLAU: A Scalable Deep Learning Accelerator Unit on FPGA

As the emerging field of machine learning, deep learning shows excellent ability in solving complex learning problems. However, the size of the networks becomes increasingly large scale due to the demands of the practical applications, which poses significant challenge to construct a high performance implementations of deep learning neural networks. In order to improve the performance as well to maintain the low power cost, in this paper we design DLAU, which is a scalable accelerator architecture for large-scale deep learning networks using FPGA as the hardware prototype. The DLAU accelerator employs three pipelined processing units to improve the throughput and utilizes tile techniques to explore locality for deep learning applications. Experimental results on the state-of-the-art Xilinx FPGA board demonstrate that the DLAU accelerator is able to achieve up to 36.1x speedup comparing to the Intel Core2 processors, with the power consumption at 234mW.

I. INTRODUCTION
I N the past few years, machine learning has become pervasive in various research fields and commercial applications, and achieved satisfactory products. The emergence of deep learning speeded up the development of machine learning and artificial intelligence. Consequently, deep learning has become a research hot spot in research organizations [1]. In general, deep learning uses a multi-layer neural network model to extract high-level features which are a combination of lowlevel abstractions to find the distributed data features, in order to solve complex problems in machine learning. Currently the most widely used neural models of deep learning are Deep Neural Networks (DNNs) [2] and Convolution Neural Networks (CNNs) [3], which have been proved to have excellent capability in solving picture recognition, voice recognition and other complex machine learning tasks.
However, with the increasing accuracy requirements and complexity for the practical applications, the size of the neural networks becomes explosively large scale, such as the Baidu Brain with 100 Billion neuronal connections, and the Google cat-recognizing system with 1 Billion neuronal connections. The explosive volume of data makes the data centers quite power consuming. In particular, the electricity consumption of data centers in U.S. are projected to increase to roughly 140 C. Wang, Q. Yu, L.Gong, X.Li and X.Zhou are with University of Science and Technology of China, Hefei, 230027, Anhui, China. E-mail: {cswang,llxx,xhzhou}@ustc.edu.cn, yuiq@mail.ustc.edu.cn).
Y. Xie is with University of California at Santa Barbara, 93106, United States, E-mail:yuanxie@ece.ucsb.edu.
billion kilowatt-hours annually by 2020 [4]. Therefore, it poses significant challenges to implement high performance deep learning networks with low power cost, especially for largescale deep learning neural network models.  to operate different sizes of tile data to leverage the trade-offs between speedup and hardware costs. Consequently the FPGA based accelerator is more scalable to accommodate different machine learning applications. 2. The DLAU accelerator is composed of three fully pipelined processing units, including TMMU, PSAU, and AFAU. Different network topologies such as CNN, DNN, or even emerging neural networks can be composed from these basic modules. Consequently the scalability of FPGA based accelerator is higher than ASIC based accelerator.

II. TILE TECHNIQUES AND HOT SPOT PROFILING
Restricted Boltzmann Machines (RBMs) have been widely used to efficiently train each layer of a deep network. Normally a deep neural network is composed of one input layer, several hidden layers and one classifier layer. The units in adjacent layers are all-to-all weighted connected. The prediction process contains feedforward computation from given input neurons to the output neurons with the current network configurations. Training process includes pre-training which locally tune the connection weights between the units in adjacent layers, and global training which globally tune the connection weights with Back Propagation process.
The large-scale deep neural networks include iterative computations which have few conditional branch operations, therefore they are suitable for parallel optimization in hardware. In this paper we first explore the hot spot using the profiler. Results in Fig. I illustrates the percentage of running time including Matrix Multiplication (MM), Activation, and Vector operations. For the representative three key operations: feed forward, Restricted Boltzmann Machine (RBM), and back propagation (BP), matrix multiplication play a significant role of the overall execution. In particular, it takes 98.6%, 98.2%, and 99.1% of the feed forward, RBM, and BP operations. In comparison, the activation function only takes 1.40%, 1.48%, and 0.42% of the three operations. Experimental results on profiling demonstrate that the design and implementation of MM accelerators is able to improve the overall speedup of the system significantly.
However, considerable memory bandwidth and computing resources are needed to support the parallel processing, consequently it poses a significant challenge to FPGA implementations compared with GPU and CPU optimization measures. In order to tackle the problem, in this paper we employ tile techniques to partition the massive input data set into tiled subsets. Each designed hardware accelerator is able to buffer the tiled subset of data for processing. In order to support the large-scale neural networks, the accelerator architecture are reused. Moreover, the data access for each tiled subset can run in parallel to the computation of the hardware accelerators.

Require:
Ni: the number of the input neurons No: the number of the output neurons Tile Size: the tile size of the input data batchsize: the batch size of the input data for n = 0; n < batchsize; n + + do ; end if end for end for end for end for In particular, for each iteration, output neurons are reused as the input neurons in next iteration. To generate the output neurons for each iteration, we need to multiply the input neurons by each column in weights matrix. As illustrated in Algorithm 1, the input data are partitioned into tiles and then multiplied by the corresponding weights. Thereafter the calculated part sum are accumulated to get the result. Besides the input/output neurons, we also divided the weight matrix into tiles corresponding to the tile size. As a consequence, the hardware cost of the accelerator only depends on the tile size, which saves significant number of hardware resources. The tiled technique is able to solve the problem by implementing large networks with limited hardware. Moreover, the pipelined hardware implementation is another advantage of FPGA technology compared to GPU architecture, which uses massive parallel SIMD architectures to improve the overall performance and throughput. According to the profiling results depicted in Table I, during the prediction process and the training process in deep learning algorithms, the common but important computational parts are matrix multiplication and activation functions, consequently in this paper we implement the specialized accelerator to speed up the matrix multiplication and activation functions. In particular, the DLAU accelerator architecture has following key features:

III. DLAU ARCHITECTURE AND EXECUTION MODEL
FIFO Buffer: Each processing unit in DLAU has an input buffer and an output buffer to receive or send the data in FIFO. These buffers are employed to prevent the data loss caused by the inconsistent throughput between each processing unit.
Tiled Techniques: Different machine learning applications may require specific neural net-work sizes. The tile technique is employed to divide the large volume of data into small tiles that can be cached on chip, therefore the accelerator can be adopted to different neural network size. Consequently the FPGA based accelerator is more scalable to accommodate different machine learning applications.
Pipeline Accelerator: We use stream-like data passing mechanism (e.g. AXI-Stream for demonstration) to transfer data between the adjacent processing units, therefore TMMU, PSAU, and AFAU can compute in streaming-like manner. Of these three computational modules, TMMU is the primary computational unit, which reads the total weights and tiled nodes data through DMA, performs the calculations, and then transfers the intermediate Part Sum results to PSAU. PSAU collects Part Sums and performs accumulation. When the accumulation is completed, results will be passed to AFAU. AFAU performs the activation function using piecewise linear interpolation methods. In the rest of this section, we will detail the implementation of these three processing units respectively.

A. TMMU architecture
Tiled Matrix Multiplication Unit (TMMU) is in charge of multiplication and accumulation operations. TMMU is specially designed to exploit the data locality of the weights and is responsible for calculating the Part Sums. TMMU employs an input FIFO buffer which receives the data transferred from DMA and an output FIFO buffer to send Part Sums to PSAU.  Fig. 2 illustrates the TMMU schematic diagram, in which we set tile size=32 as an example. TMMU firstly reads the weight matrix data from input buffer into different BRAMs in 32 by the row number of the weight matrix (n=i%32where n refers to the number of BRAM, and i is the row number of weight matrix). Then, TMMU begins to buffer the tiled node data. In the first time, TMMU reads the tiled 32 values to registers Reg a and starts execution. In parallel to the computation at every cycle, TMMU reads the next node from input buffer and saves to the registers Reg b. Consequently the registers Reg a and Reg b can be used alternately.
For the calculation, we use pipelined binary adder tree structure to optimize the performance. As depicted in Fig.  2, the weight data and the node data are saved in BRAMs and registers. The pipeline takes advantage of time-sharing the coarse-grained accelerators. As a consequence, this implementation enables the TMMU unit to produce a Part Sum result every clock cycle.

B. PSAU architecture
Part Sum Accumulation Unit (PSAU) is responsible for the accumulation operation. Fig. 3 presents the PSAU architecture, which accumulates the part sum produced by TMMU. If the Part Sum is the final result, PSAU will write the value to output buffer and send results to AFAU in a pipeline manner. PSAU can accumulate one Part Sum every clock cycle, therefore the throughput of PSAU accumulation matches the generation of the Part Sum in TMMU.

C. AFAU architecture
Finally, Activation Function Acceleration Unit (AFAU) implements the activation function using piecewise linear interpolation (y=ai*x+bi, x∈[x 1 ,x i+1 )). This method has been widely applied to implement activation functions with negligible accuracy loss when the interval between x i and x i+1 is insignificant. Eq. (1) shows the implementation of sigmoid function. For x>8 and x≤-8, the results are sufficiently close to the bounds of 1 and 0, respectively. For the cases in -8<x≤0 and 0<x≤8, different functions are configured. In total we divide the sigmoid function into four segments.
Similar to PSAU, AFAU also has both input buffer and output buffer to maintain the throughput with other processing units. In particular, we use two separate BRAMs to store the values of a and b. The computation of AFAU is pipelined to operate sigmoid function every clock cycle. As a consequence, all the three processing units are fully pipelined to ensure the peak throughput of the DLAU accelerator architecture.

IV. EXPERIMENTS AND DATA ANALYSIS
In order to evaluate the performance and cost of the DLAU accelerator, we have implemented the hardware prototype on the Xilinx Zynq Zedboard development board, which equips ARM Cortex-A9 processors clocked at 667MHz and programmable fabrics. For benchmarks, we use the Mnist data set to train the 784×M×N×10 Deep Neural Networks in Matlab, and use M×N layers weights and nodes value for the input data of DLAU. For comparison, we use Intel Core2 processor clocked at 2.3GHz as the baseline.
In the experiment we use Tile size=32 considering the hardware resources integrated in the Zedboard development board. The DLAU computes 32 hardware neurons with 32 weights every cycle. The clock of DLAU is 200MHz (one cycle takes 5ns). Three network sizes-64×64, 128×128, and 256×256 are tested.

A. Speedup Analysis
We present the speedup of DLAU and some other similar implementations of the deep learning algorithms in Table  II. Experimental results demonstrate that the DLAU is able to achieve up to 36.1x speedup at 256×256 network size. In comparison, Ly&Chows work [5] and Kim et.als work [7] present the work only on Restricted Boltzmann Machine algorithms, while the DLAU is much more scalable and flexible. DianNao [6] reaches up to 117.87x speedup due to its high working frequency at 0.98GHz. Moreover, as DianNao is hardwired instead of implemented on a FPGA platform, therefore it cannot efficiently adapt to different neural network sizes. Fig. 4 illustrates the speedup of DLAU at different network sizes-64×64, 128×128, and 256×256 respectively. Experimental results demonstrate a reasonable ascendant speedup with the growth of neural networks sizes. In particular, the speedup increases from 19.2x in 64×64 network size to 36.1x at the 256×256 network size. The right part of Fig. 4 illustrates  how the tile size has an impact on the performance of the DLAU. It can be acknowledged that bigger tile size means more number of neurons to be computed concurrently. At the network size of 128×128, the speedup is 9.2x when the tile size is 8. When the tile size increases to 32, the speedup reaches 30.5x. Experimental results demonstrate that the DLAU framework is configurable and scalable with different tile sizes. The speedup can be leveraged with hardware cost to achieve satisfying trade-offs.

B. Resource utilization and Power
Table III summarizes the resource utilization of DLAU in 32×32 tile size including the BRAM resources, DSPs, FFs, and LUTs. TMMU is much more complex than the rest two hardware modules therefore it consumes most hardware resources. Taking the limited number of hardware logic resources provided by Xilinx XC7Z020 FPGA chip, the overall utilization is reasonable. The DLAU utilizes 167 DSP blocks due to the use of the Floating-point addition and the Floatingpoint multiplication operations.   [7], the BRAM utilization of DLAU is insignificant. This is due to the tile techniques so that large scale neural networks can be divided into small tiles, therefore the scalability and flexibility of the architecture is significantly improved. In order to evaluate the power consumption of accelerator, we use Xilinx Vivado tool set to achieve power cost of each processing unit in DLAU and the DMA module. The results in Table IV-B depict that the total power of DLAU is only 234mW, which is much lower than that of DianNao (485mW). The results demonstrate that the DLAU is quite energy efficient as well as highly scalable compared to other accelerating techniques. To compare the energy and power between FPGA based accelerator and GPU based accelerators, we also implement a prototype using the state-of-theart NVIDIA Tesla K40c as the baseline. K40c has 2880 stream cores working at peak frequency 875MHz, and the Max Memory Bandwidth is 288 (GB/sec). In comparison, we only employ 1 DLAU on the FPGA board working at 100MHz. In order to evaluate the speedup of the accelerators in a real deep learning applications, we use DNN to model 3 benchmarks, including Caltech101, Cifar-10, and MNIST, respectively. Fig. 5 illustrates the comparison between FPGA based GPU+cuBLAS implementations. It reveals that the power consumption of GPU based accelerator is 364 times higher than FPGA based accelerators. Regarding the total energy consumption, the FPGA based accelerator is 10x more energy efficient than GPU, and 4.2x than GPU+cuBLAS optimizations.
Finally Fig. 6 illustrates the floor plan of the FPGA chip. The left corner depicts the ARM processor which is hardwired in the FPGA chip. Other modules, including different components of the DLAU accelerator, the DMA, and memory interconnect, are presented in different colors. Regarding the programming logic devices, TMMU takes most of the areas as it utilizes a significant number of LUTs and FFs.

V. CONCLUSION AND FUTURE WORK
In this article we have presented DLAU, which is a scalable and flexible deep learning accelerator based on FPGA. The DLAU includes three pipelined processing units, which can be reused for large scale neural networks. DLAU uses tile techniques to partition the input node data into smaller sets and compute repeatedly by time-sharing the arithmetic logic. Experimental results on Xilinx FPGA prototype show that DLAU can achieve 36.1x speedup with reasonable hardware cost and low power utilization.
The results are promising but there are still some future directions, including optimization of the weight matrix and memory access. Also the trade-off analysis between FPGA and GPU accelerators is another promising direction for large scale neural networks accelerations.


Exploiting Convolution Filter Patterns for Transfer Learning

In this paper, we introduce a new regularization technique for transfer learning. The aim of the proposed approach is to capture statistical relationships among convolution filters learned from a well-trained network and transfer this knowledge to another network. Since convolution filters of the prevalent deep Convolutional Neural Network (CNN) models share a number of similar patterns, in order to speed up the learning procedure, we capture such correlations by Gaussian Mixture Models (GMMs) and transfer them using a regularization term. We have conducted extensive experiments on the CIFAR10, Places2, and CMPlaces datasets to assess generalizability, task transferability, and cross-model transferability of the proposed approach, respectively. The experimental results show that the feature representations have efficiently been learned and transferred through the proposed statistical regularization scheme. Moreover, our method is an architecture independent approach, which is applicable for a variety of CNN architectures.

Introduction
The CNN models are found to be successful at various computer vision tasks such as image classification [14,27,28], object detection [22], image segmentation [19,21], and face recognition [20], where large-scale datasets [24,32,20] are available. Nevertheless, the performance of CNN models significantly reduces, when training data is limited or the domain of the training set is far from the test set. Today, the most successful and practical solution to address lack of annotated large dataset is training the networks on large-scale annotated datasets like ImageNet [24] and Places [32], then finetuning these pre-trained networks for specific problems. Thanks to community, the pretrained models of well-known architectures like AlexNet [18], VGG-16 [27], GoogLeNet [28], and ResNet [14] can be found available online. However, when some architectural changes are needed, these pre-trained networks can-not be used. For such cases, it is necessary to train models on large datasets, then, finetune for the particular problem. Unfortunately, while getting nearly human performances on a lot of applications with these models, training these networks on large datasets is still a significant problem and a very time consuming process.
With recent advances in deep learning such as Residual Learning [14], successful networks have become more and more deep, and training these models become harder in terms of complexity and time. To find a solution for this problem, inspired by What makes a good detector? [11], we have investigated What makes a good CNN filter? Two successful CNN models, AlexNet [18] and VGG-16 [27], are analyzed from a statistical perspective, and we realized that these models show similar patterns and redundancies. In addition to our findings, in [7], authors show that 95% of weights of neural networks could be predicted without any reduction in accuracy. This leads us to the idea that we can use these redundancies and patterns for learning better representations quickly by transfer learning. Similar to the methods used in [2,4,11], we introduce a regularization term for transferring the statistical information to improve the learning scheme. First, the statistical distribution of convolution filters from a well-trained network is learned with a Gaussian Mixture Model. Next, the newly trained model is encouraged to show similar statistics with source models using the regularization term. Extensive experiments on the CIFAR10 [17], Places2 [32], and CMPlaces [4] datasets show that the proposed approach is generalizable and the networks can quickly learn a representation with statistical regularization, which could efficiently be transferred to another task and cross-domains. The overview of our proposed method can be seen in Figure 1. The rest of the paper is organized as follows: in Section 2 related works are summarized, in Section 3 our detailed statistical analysis is reviewed and the regularization term is introduced, in Section 4 experimental results are presented and discussed, and finally in Section 5 the paper is concluded. Figure 1: Overview of the proposed system. The blue (top) network is a well-trained CNN where its distribution is learned via GMM. When new red (bottom) network being trained, along with the classification loss, the new regularization term λ · R(w) that measures the statistical difference of two distributions are minimized. The statistical knowledge could be transferred from source CNN to target CNN with our regularization term. Best viewed in color.

Related Work
Network Distillation: The aim of network distillation approaches are transforming larger networks to smaller ones, while not losing so much information that was learned in the large network. The pioneering work is conducted by Bucila et al. [3], where their aim is to compress the ensemble of models to a single model without significant accuracy loss. Later, Hinton et al. [15] optimized the smaller network to show similar softmax output of cumbersome model. Then, Romero et al. [23] suggest that in addition to softmax outputs, intermediate representations could be used for distilling the network. Recently distillation is also applied in the cross-modal settings [13,1]. The major drawback of these methods is the necessity to train a large network before using it to train a smaller one. Also these models mainly match the outputs of the networks, whereas we regularize the internal weights of the network.
Domain Adaptation: Domain adaptation is the problem of learning a model that generalizes target domain examples besides source domain ones while learning only from source domain. With the success of CNNs, domain adaptation works have been focused on CNNs and several successful methods have been proposed. For instance, Ganin & Lempitsky [9] introduced gradient reversal layer which act like identity transform in forward pass of CNN and change the sign of the gradient and scale in the backward pass. In their work, they added a domain classifier to the end of the feature map and tried to predict the domain of examples. In backpropagation, they changed the gradient using gradient reversal layer and forced the networks to learn domain invariant features by maximizing loss of domain classifier. Tzeng et al. [29] added new terms into objective function of CNN to both increasing domain confusion and transferring inter-class knowledge. The first term domain confusion loss forces the networks to learn domain invariant features and soft label loss forces the feature of the same class to be similar for both source and target domain. In contradistinction to Ganin & Lempitsky [9], some target labels must be available for optimizing soft label loss. Moreover, recently some other works [10,12,30,26] focused on domain adaptation problems for CNNs. While these methods focus on transferring information about structure of data, our method focuses on transferring more local information.
Statistical Transfer: Statistical transfer is learning statistical properties from a source and to use this statistical knowledge for improving learning procedure. For instance, Aytar & Zisserman [2] proposed part-level transfer regularization which transfers parts of source detectors instead of whole detector. Additionally, they take advantage of part co-occurrence statistics. For example, if there are wheels in the picture, probably another wheel would also appear. They calculated these co-occurrence statistics using the source data and transfer these statistics when a new object detector is learned. Moreover, in the era of hand crafted features, Gao et al. [11] analyzed famous HOG [6] templates of successful object detectors and made two observations. Firstly, they observed that activations of individual cell models had some correlations, secondly local neighborhoods of cells also showed the same characteristic. Furthermore, since they wanted to transfer local information in contrast to templates like in [2], they defined their priors such that the correlations from the source model could be transferred to target model without global template align-ment. In a recent work [4], GMMs are used for aligning cross-modal data. In this approach, statistics of activation maps of different layers for a modality are learned, and the other modalities are forced to show similar statistics in the activations via a regularization term. The proposed method is capable of aligning modalities by using this regularization term when there is no strong alignment between modalities.
Our work is influenced by Gao et al. [11] and Castrejon et al. [4]. The first work directs us to analyze the weights of convolution filters and find the correlation between the filters, and the other one to use the GMM to enforce them to show similar statistics in a non-convex optimization problem. Different from these works, our regularization forces the weights to show similar statistics and transfers local correlations of convolution filters.

Statistical Regularization
In this section, first, we present our analysis on CNN models from a statistical perspective. Next, we describe our approach for capturing statistical knowledge from a CNN and present how to transfer this knowledge via a regularization term.

Statistical Analysis
For investigating the statistical behavior of a CNN model, we have used VGG-16 [27] and AlexNet [18] models that were trained on ImageNet [24] and Places-365 [32] datasets. Especially, weights of convolution filters are analyzed for the four models. We tried to answer following three questions: (i) are filters separable into clusters?; (ii) how much similar are the filters inside a cluster?; and (iii) how all filters are distributed over clusters?. In order to obtain this information, all 3 × 3 filters of a model are clustered to ten different clusters using the k-means algorithm and each cluster's covariance matrix and mean value were calculated and visualized. The covariance matrices would provide information about how members of a cluster are correlated. Mean values show whether the sets are similar to each other or not. For instance, it can be seen in Figure 2 (a) the covariance matrices of all the four models show that there is some shared behavior of learned filters across layers and architectures. Also, the cluster centers depict different similarities. Generally, all models have clusters that their mean values are accumulated on the left, right, top, and bottom. Moreover, when we look at the distribution of filters to clusters, the models show different characteristics. While the distributions are Gaussian-like in VGG-16 models, the distributions cannot be fitted to a known distribution in the AlexNet models and most of the clusters have roughly the same number of filters. However, both AlexNet and VGG-16 models show similar distributions by model wise while they are trained on different dataset. The mean values are shown in Figure 2 (b), and the distributions can be seen in Figure 3.

The Proposed Method
Since our aim is capturing statistical properties of "good" convolution filters and transfering this statistical knowledge to another network, we modeled the convolution filters' distributions and enforced the network to show similar distributions. Similar to our work, [4] forces the activations of networks to show similar distributions across modalities for aligning cross modal data. In [4], the authors used both mixture and single Gaussian distributions for modelling activations, and the mixture models outperformed single Gaussian model in their problem. Since, our aim is also similar to their work -capturing statistical knowledge and transfering it-we have decided to use mixture models for modeling distributions. The main difference is that we force the weights of filters instead of the activations to show similar distributions across the networks. Let x n and y n be a training image and its corresponding label. We want to minimize where z n (x n ,w) is output of the network. We have added a regularization term R to the loss term that represents the negative log likelihood of a convolution filter to encourage the network to learn the weights that are statistically similar to another network. For filter w i and distribution P we define R such that, The distribution P is modeled as GMM, therefore P would be, where K is the number of mixtures and k π k = 1 , π k ≥ 0 ∀ k . Therefore, the total negative log likelihood for N convolution filters can be defined as, where N is specified as, The derivate of the negative log likelihood must be calculated exactly or analytically, since the network is trained with back-propagation. Still, calculating exact derivative in every iteration would be very expensive during training,  therefore, we approximately calculate the derivative. To approximately calculate the derivative of a convolution filter, we first pick the mixture component N (µ s , Σ s ) that the probability of the convolution filter is maximized. Next, the derivate is calculated using that single Gaussian. The partial derivative of R with respect to a convolution filter w i would be Finally, our complete loss term is defined as where the first and second terms are classification loss and weight decay terms, and the last one is our regularization term. The λ is a hyperparameter that controls the regularization. For experiments in this paper, we have used distribution P , which is learned from convolution filters of VGG-16 model trained on ImageNet using Expectation maximization (EM) algorithm with K = 1000 components. To reduce the number of parameters, we assumed that the covariances Σ k are diagonal. Also, K-Means algorithm was run on filters for decreasing the convergence time of EM algorithm. While it is known that, different layers, especially first layers, show different characteristics than other layers, we have used a single GMM -P -for all filters. Because, if multiple distributions are used, a layer alignment relationship between layers are needed in transfer time.

Experiments and Results
Our hypothesis was that well-trained CNN models show similar statistical patterns, and we could exploit this information in the training phase. Our experiments show that we can learn a better representation faster with the statistical regularization. For example, in section 4.1, we validate that with regularization the convolution filters create more general representations and show less overfitting characteristics than regular trained filter's representations. Also in section 4.2, we see how the regularization helps to learn representations that are successfully transferred to another task. Furthermore, while transferring across task and domains is widely studied, transfer learning for cross modal data is not a well studied problem. We show that our method could be applied to cross-modal data as well. In section 4.3, we present our experiments applied to a cross-modal dataset CMPlaces [4].

Generalization
In this section, we evaluate if the statistical regularization helps the generalizability of the learned representation. For this purpose, we have used middle sized CIFAR-10 dataset [17] and CNN architecture described in [25]. Firstly, we trained the network with and without regularization and stopped the training at 10k, 15k, 20k, and 25k iterations. Later on we freeze all 3 × 3 filters and continue learning with training data. Since only the last layers change during the training and the features extracted from convolution filters are not so generic, the validation loss starts to increase and the validation accuracy is dropped.
We compared the networks, whose training stopped in different iterations with and without regularization. In our experiments, the training loss is started to oscillate in a small interval and does not change much, that is, neither increases nor decreases. However, the validation loss changed. The gap between training and test loss indicates how generalizable the network is.
When we compared the models initialized at 10k, 15k, 20k, and 25k iterations, always the test loss of regularized versions increased slowly compared to those of nonregularized versions. For instance, as can be seen in Figure  4 (a,b,c,d), the difference between training and test loss increases slowly in regularized models than non-regularized models. Interestingly, as regular training time (no freezing) increases, the gap between the regularized and normal training reduces.

Task Transfer
In this section, we show that the quickly learned representations can also be transferred to another task more successfully with regularization. To validate our claim, we try to transfer filter distributions from ImageNet to Places2 [32] dataset. The trivial solution for a classification problem is training the network on ImageNet and finetuning the models on the new task afterwards. We follow the same procedure in this section for our experiments, but we want to finish the pre-training stage as early as possible. Also we want to show that our regularization can be used with finetuning. To evaluate our method's performance, we first train VGG-F model introduced in [5] on ImageNet [24] data with and without regularization. As in Section 4.1, we take snapshots from 10k, 25k, and 50k iterations. Next, we start finetuning on Places2 data [32]. When we compare regularization effect, we see that the regularization can help to learn better representations in the early iterations. For example, as can be seen in Figure 5 (a), the performance of the models that are initialized at the weights learned in pre-training with only 10k iterations, the regularized version shows a better performance than the non-regularized one. When we examine the models initialized at 25k iterations, the performance difference between regularized and non-regularized versions reduces. Finally, in 50k iterations, nearly there is no difference in the performance. The test loss/iterations plots are shown in Figure 5 (a), (b), (c). This experiment shows that as pre-training time increases, the gain obtained from regularization decreases. However, for limited amount of pre-training time, the regularization could increase the efficiency of the pre-training.

Modality Transfer
While CNNs performances are very good at various computer vision tasks for real-world images, most of the computer vision algorithms fail at non-real images. This shows that generalization performance of the computer vision algorithms is not good for cross modal data. Some recent works [4,8] focused on this problem and recently Castrejon et al. [4] has introduced a new Cross-Modal dataset. In the dataset there are five different modalities for each scene type such as natural image, sketch, clip art, spatial text, and description. As in section 4.2, we firstly train VGG-F models on a large dataset -Places2-and finetune on CMPlaces data. We have used clipart and sketch data to  evaluate our performance, since we are interested in the visual domain. We take snapshots from 25k and 50k iterations from the regularized and non-regularized networks. Next, we start to finetune on sketch and clipart data and compare their accuracies. Also, we finetune VGG-F model converged on ImageNet data and compare with our pre-trained models. Similar to our experiments, by rising the training pre-training time the performance increases. Also, the regularization helps to learn better representations and increases the performance of pre-training for finetuning. Although, there is a significant gap between the converged ImageNet model and our regularization in the sketch data, there is not a substantial difference between 50k iteration with regularization and the converged ImageNet model on clipart modality. The results for both modalities can be seen in Table 1 and Table 2. These results show that instead of training ImageNet until the model converges, we could train the models using regularization with only few iterations and could employ these pre-trained networks for cross-modality transfer.

Implementation Details
We have used Caffe [16] deep learning framework in our experiments. Moreover, we have implemented our special convolution layer for applying statistical regularization.  When VGG-F model is trained on ImageNet and Places datasets, stochastic gradient descent with 0.01 learning rate is used for optimization. In the CIFAR experiments, we have used the same parameters described in [25]. Finally, the Gaussian mixture models are learned using the VLFeat library [31].

Conclusion
In this paper, we analyzed convolution filters of wellknown CNN architectures and found that they share a number of common patterns and redundancies that could be exploited for transfer learning. Gaussian Mixture Models are used for capturing these statistical patterns and a new regularization term is introduced for transferring such patterns to other networks. Our experiments show that we could learn good representations that are transferable to the other tasks and cross-domains quickly with regularization. For instance, we achieved around 25% improvement on the sketch modality in the cross-modal dataset under limited pre-training time. Also, our method gets similar performance on clipart data with converged model that pre-trained on the ImageNet, while pre-training stopped at 50k iterations in our method.


Modeling an impact of road geometric design on vehicle energy consumption

Some roads connect traffic origins and destinations directly, some use winding, indirect routes. Indirect connections result in longer distances driven and increased fuel consumption. A similar effect is observed on congested roads and mountain roads with many changes in altitude. Therefore a framework to assess road networks based on energy consumption is proposed. It has been shown that road geometry has significant impact on overall traffic energy consumption and emissions. The methodology presented in the paper analyzes impact of traffic volume, shares of vehicle classes, road network configuration on the energy used by the vehicles. It can be used to optimize energy consumption with efficient traffic management and to choose optimum new road in the design phase. This is especially important as the energy consumed by the vehicles shortly after construction supersedes the energy spent for the road construction.

Introduction
The European Union puts a lot of effort on policies for energy efficiency and to decrease emissions in transport [1]. Transport sector is contributing nearly a third of CO2 emissions and energy consumption within the EU [2]. Most efforts to combat this problem are directed towards traffic mode shift, use of energy-efficient vehicles and alternative fuels [3]. As a result there was a decline of overall emissions by 3% in 2012 and by 2014 average light vehicle emissions were below the targets set for 2015. These values were obtained in part as the EU requires member states to record data on registering motor vehicles [4]. Similar policies are being practiced elsewhere in the world, such as in the USA and China [5,6]. Such transport policies also focus on limiting road traffic in cities, increasing use of public transport and on shifting from road to other modes of transport. Such approaches work well in densely populated areas, but in sparsely populated areas road traffic is still, necessarily, the transportation mode of choice. An example of a study of reducing emissions in a densely populated area is the article [7], where a case study from London is described. An alternative to such approaches can be better traffic management that could result in less congestion and lower emissions.
When roads are preferred option and shift to other modes is difficult to achieve, more focus should be put on relation between greenhouse emissions and road network configuration including traffic management. Traffic simulation and emissions models are a tools that should be applied to optimize the traffic network.
Traffic microsimulation models need to be coupled with emission models in order to obtain individual and aggregated vehicle emissions. Vehicle emission models can be divided into microscopic and macroscopic. Macroscopic emission models are used to estimate emissions based on traffic volume and average speed, ignoring instantaneous accelerations. More detailed are the microscopic emission models that rely on inputs such as speed, acceleration, vehicle weight and road grade.
Regardless if the model is microscopic or macroscopic, most models provide no direct physical model of fuel consumption, rather tables or functions that are based on statistical analysis of different variables that have an impact on emissions [8,9]. One example of a macroscopic model is the ARTEMIS study that involves CO2 and NOX emissions dependant on vehicle speed and driving cycles [10,11]. A statistical analysis of the dependence of truck fuel consumption on speed and road grade was conducted in an experimental study by [12] conducted on highways in realistic traffic conditions. Fuel volume flow was measured using the vehicles' CAN (Controller Area Network) bus interface. Such results can be used for comparative purposes and for model tuning but not for direct estimation of instantaneous fuel consumption for different types of vehicles that is required for microsimulation. For general energy demand forecasting models with lower level of detail are used. In [13], where a long-term forecasting was simulated, a straight-forward model based on mileage travelled and average fuel economy was applied.
Notable microsimulation models are Versit+, PHEM, CMEM and VT-Micro. Versit+ is a model that is based on a vehicle type, velocity and acceleration [14] that was developed by the TNO from Netherlands. It has interface with AIMSUN traffic microsimulation software. The AIMSUN is based on the Gipps model therefore velocity must be differentiated to obtain acceleration [15]. PHEM is an emissions model that was developed by the University of Graz and has interface for the SUMO microscopic traffic simulator [16][17][18]. CMEM and VT are models that were development in the USA. The CMEM was developed during the 90s and was afterwards, while the VT (Virginia Tech) is a newer model. Both of them are based on velocity, acceleration and factors that depend on vehicle types.
Since new fuels and electric vehicles are becoming increasingly popular, calculating energy impact of traffic can be difficult. A lot more straight-forward approach would be to rely on the energy required by the vehicle directly. Among all approaches, the VSP (vehicle specific power) model was the only one that was directly based on power required by the vehicle [19,20]. Vehicle specific power is based on energy conversions, changes of potential and kinetic energy. Power needed for acceleration results in a change of kinetic energy and is therefore affected by congestion; road grade, too, has an impact on potential energy change. Apart from the energies, air and rolling resistances are also modelled. There are two frequently mentioned formulas for VSP in the literature. The first was developed by Jimenez and is based directly on changes of energy, and further derived to include drag and rolling resistance coefficients and accelerations [19]. Another formula was developed by Zhai; this also introduces several coefficients [20]. The main difference between the two is that Jimenez uses wind velocity as an independent variable, which is the reason this formula is preferred; in case of unidirectional tunnels the air velocity is usually in the direction of travel, which results in lower air resistance as opposed to the open road.
The downside of the VSP model is that it does not account for engine idling, which can be significant in case of congested traffic. According to [21] and [22], fuel consumption during engine idling is between 0.5 and 1.5 l/h for light vehicles and a bit over 3 l/h for heavy vehicles. The VSP is also suitable for adaptation to analysis of alternative drivetrains as it is based on energy consumption with a clear physical background. It is not a regressive model of fuel use. It merely estimates power at the wheels. The power at the engine should be higher 2 E3S Web of Conferences 22, 00107 (2017) DOI: 10.1051/e3sconf/20172200107 ASEE17 as transmission and engine efficiencies should be considered as well. A relationship between the vehicle specific power and fuel consumption rate (flow) was conducted by [23,24].
An exhaustive study involving vehicle efficiency defines separate ratios for powertrain, vehicle and primary energy efficiency [25]. It can be deduced that due to many different fuels that are being used, it is more feasible to have a model that is based on energy and then to calculate primary energy use from the energy at the wheels.

Extraurban case study -Karawanks tunnel
An analysis of the Karawanks road tunnel connecting northern Slovenia with southern Austria was carried out. It is a single tube tunnel with bidirectional traffic connecting motorways on both sides of the border. It is located under the Karawanks mountain ridge, which is approximately 2500 m high. This particular route is interesting as all alternative routes are much longer and include many serpentine segments and changes in elevation (as can be seen in figures 3 and 4). In figure 4 the axis represents the geographical longitude and latitude. The road coordinates were obtained from OpenStreetMap, while the elevation data are from GMTED 2010 database [26]. It is evident that the tunnel route (black) is a lot shorter than the alternative route because it crosses the mountain ridge (the white route in figure 4). It is, however, still longer than the hypothetical ideal route (dashed line). Length and elevation comparison of the routes are shown in table 1.   It would be possible to split traffic data into several time periods, but since detailed traffic data were not available analysis using an average number of vehicles per day was conducted: separating them into personal vehicles, busses, light commercial vehicles and heavy goods vehicles. Vehicle weight for each class was approximated. An overview of traffic data can be seen in the table 2. As can be seen, the traffic volume is relatively low with a high share of personal vehicles.
After obtaining the traffic data, we assigned each vehicle class an approximate weight and length. The vehicle length is needed for microsimulation as well. After obtaining vehicle trajectories, VSP calculation along the route was carried out. The VSP is highly dependent on road grade and accelerations. Changes in road grade have many peaks that are not evident in reality because of the resolution of satellite elevation data and the very steep changes in elevation in these narrow Alpine valleys. These peaks will not necessarily yield wrong results as there are significant changes in elevation on the actual road which are spread over longer distances compared to those of the satellite based elevation data. If the analysed road was on a flat area there such irregularities would be minimized.
Calculated average values for energy consumption that were obtained by pre-processing microsimulation vehicle trajectories of each vehicle in a class are shown in table 4, where route 1 is the tunnel route, 2 is the alternative route and 3 is the ideal route. Based on data from table 4 the daily amount of energy consumed by all vehicles on all 3 routes can be calculated according to the eq. 7. The results can be seen in table 6, which shows that the tunnel route is about 3 times more efficient than the alternative route (230% more energy is consumed). An interesting overview is the average daily power needed to drive all the vehicles that pass on the route for all three options. The difference between the tunnel and alternative routes is 18 MW which is an order of magnitude of a smaller electrical power station. This means that the impact of road geometry on fuel consumption is far from negligible. I

Urban case study -Celovska street, Ljubljana
Celovska street is one of major streets in the city of Ljubljana that connects northern bypass ring motorway with city centre. Most of the street has 2 lanes per direction and includes several semaphored intersections. Hourly traffic volume values were obtained from city authorities and are based on inductive loop counters on the street [27]. Average monthly hourly traffic volume from bypass towards the centre can be seen in figures 4 and 5.
A simulation was carried out using predefined flow of moderate 400 vehicles per hour which is normal daily traffic (excluding peaks) plus random trips across the network. 3 different traffic light control programmes were evaluated among which are: 1fixed 90 second cycle programme for all traffic lights, 2adaptive programme with extensions of priority phases and maximum allowed cycle duration of 125 s and 2an adaptive programme with phase extensions combined with forced phase switching and synchronisation among junctions. Goal of this simulation is to assess impact of energy consumption on different traffic light management programmes. It is evident that according to the model, there are dramatic differences between noncongested traffic (programmes 2 and 3) and congested (programme 1).  control programmes. These results are still preliminary as more detailed traffic data are needed. It should be noted that modelled traffic control programmes are not currently used programmes but are an evaluation of different possible traffic control programmes on a realistic traffic network with real traffic volume.

Conclusions
A methodology to evaluate an impact of road network configuration on energy consumed by the vehicles using the vehicle specific power model presented. It has been shown that it is possible to save enormous amounts of energy by making road networks more energy efficient, by building more direct, shorter routes without significant changes in road grade or by controlling traffic lights more efficiently. First is especially evident in mountainous areas where open roads must pass through significant variations in elevation. This results in longer routes with higher road grades, causing higher energy consumption. In the cities same effect was observed by avoiding congestions with more efficient control of traffic lights.
Vehicle energy consumption was modelled with an assumption that no braking energy is recovered. In order to simulate electric and plug-in hybrid vehicles, simple modifications to the VSP model are needed.