extfieldsofstudy:
Physics
Computer Science
Mathematics

text:
Private and Quantum Capacities of More Capable and Less Noisy Quantum Channels

Two new classes of quantum channels, which we call more capable and less noisy, are introduced. The more capable class consists of channels such that the quantum capacities of the complementary channels to the environments are zero. The less noisy class consists of channels such that the private capacities of the complementary channels to the environment are zero. For the more capable class, it is clarified that the private capacity and quantum capacity coincide. For the less noisy class, it is clarified that the private capacity and quantum capacity can be single letter characterized.

I. INTRODUCTION
One of the most important problem in quantum information theory is to determine the quantum capacity of a noisy quantum channel. The capacity is defined as the transmission rate optimized over all possible quantum error correcting codes such that decoding errors vanish in the limit of asymptotically many uses of the channel.
Mathematically, a quantum channel can be described by the trace preserving completely positive (TPCP) map from the input system to the output system. By using the Stinespring dilation of the TPCP map, we can naturally define a complementary channel to an environment system, and we can regard the noisy quantum channel as a wire-tap channel [1,2] from the sender to the legitimate receiver and the eavesdropper who can observe the environment system of the channel (eg. see [3]). Then we can define the private capacity of the noisy quantum channel as the transmission rate optimized over all possible wire-tap codes such that decoding errors and information leakage vanish in the limit of asymptotically many uses of the channel.
The private capacity and quantum capacity of noisy quantum channels were established in [4][5][6][7]. However unlike the capacity formula of a classical noisy channel or the private capacity formula of a classical wire-tap channel, the private capacity and quantum capacity formulae are not single letter characterized, i.e., the formulae involve the limit with respect to the number of channel uses, and they are not computable. Indeed, some numerical evidences clarified that the expressions in the capacity formulae are not additive [8][9][10][11], and the single letter characterization is not possible in general at least by using the same expressions.
A quantum channel is called degradable if there exists another degrading channel such that the conjunction of the channel to the legitimate receiver and the degrading channel coincide with the complementary channel to the eavesdropper. In such a case, the single letter charac- * Electronic address: shun-wata@is.tokushima-u.ac.jp terizations of the private capacity and quantum capacity were established [3,12].
A quantum channel is called conjugate degradable if there exists another degrading channel such that the conjunction of the channel to the legitimate receiver and the degrading channel coincide with the complementary channel to the eavesdropper up to complex conjugation. In such a case, the single letter characterizations were also established [13].
To date, all quantum channel whose capacities are single letter characterized are degradable or conjugate degradable 1 , and it is important to clarify a broader class of quantum channels such that the single letter characterizations are possible.
Aside from the possibility of the single letter characterizations, there is also another interesting problem. In the quantum information theory, the private information transmission and the quantum information transmission are closely related [4,[14][15][16], and the possibility of the latter implies the possibility of the former. However, the private information transmission and the quantum information transmission are not exactly equivalent. Indeed, although the private capacity and quantum capacity coincide for degradable quantum channels [17], the former can be strictly larger than the latter in general. Especially the private capacity can be positive even if the quantum capacity is zero [18]. Thus it is important to clarify a condition on quantum channels such that the private capacity and quantum capacity coincide or not.
To shed light on the above mentioned two problems, we introduce two new classes of quantum channels, which we call more capable and less noisy. The more capable class consists of channels such that the quantum capacities of the complementary channels are zero. The less noisy class consists of channels such that the private capacities of the complementary channels are zero. Later, these definitions turn out to be natural analogies of the partial orderings, more capable and less noisy, between classical channels [19].
The inclusive relation of the degradable, the conjugate degradable, the less noisy, and the more capable classes are summarized in Fig. 1. In this paper, we show that the private capacity and quantum capacity coincide for channels in the more capable class. Furthermore, we also show that the private capacity and quantum capacity can be single letter characterized for channels in the less noisy class. These results provide partial solutions to the above mentioned two problems.
The rest of the paper is organized as follows. In Section II, we review some known results on the private capacity and quantum capacity of quantum channels. In Section III, the more capable and less noisy classes are introduced, and we state our main results. In Section IV, we summarize certain properties implied by more capable and less noisy, and show proofs of our main results. We finalize the paper with conclusion in Section IV.

II. PRELIMINARIES
Let N B be a quantum channel from an input system H A to an output system H B . By using the Stinespring dilation (eg. see [3]), there exist an environment system H E and an isometry U BE from H A to the joint system for every input ρ, where Tr E is the partial trace with respect to the environment system. By using this representation, we can naturally define another channel which is usually called the complementary channel of N B . Although the Stinespring dilation is not unique, the following arguments do not depend on the choice of the dilation because two dilations can be converted to each other by applying a local unitary to the environment systems.
Throughout the paper, we basically follow the notations from [3,20]. The von Neumann entropy of a density matrix ρ is defined by H(ρ) = −Trρ log ρ, and the quantum relative entropy between ρ and σ is defined by D(ρ σ) = Trρ(log ρ − log σ). For input state ρ A to the channel N B , the coherent information is de- When the input state is clear from the context, we omit the subscript and denote I c (A B). The quantum mutual information of ρ XB on the joint system is defined by then the quantum mutual information can be written as When the legitimate receiver can observe the output of N B and the eavesdropper can observe the output of N E , the private capacity is characterized by [4,5] C p (N B ) = lim where where {ρ u A } are states (not necessarily pure states) on H A indexed by u ∈ U, and P U is a probability distribution on a finite set U.
On the other hand, when the sender want to transmit quantum information to the receiver through the channel N B , the quantum capacity is characterized by [4,6] where Both Eqs. (1) and (2) cannot be single letter characterized, i.e., C p (N B ) > C (1) p (N ) and Q(N ) > Q (1) (N ) in general [8][9][10][11]. Furthermore, the private capacity can be strictly larger than the quantum capacity, i.e., C p (N B ) > Q(N B ) for some channels [18].
The channel N B is said to be degradable if there exists a TPCP map D such that D • N B = N E . This is a quantum analogue of degraded broadcast channel [21]. When N B is degradable, then it is known [3,12] that the single letter formulae hold, i.e., C p (N B ) = C (1) Let C denote entry-wise complex conjugation with respect to a fixed basis of H E . Then, the channel N B is said to be conjugate degradable if there exists a TPCP map D such that D • N B = C • N E . When N B is conjugate degradable, it is known that Q(N B ) = Q (1) (N B ) [13]. Later, it will turn out that C p (N B ) = Q(N B ) = Q (1) (N B ) also holds.

III. MAIN STATEMENTS
In this section, we introduce two new classes of quantum channels, and show our main results.

Definition 1
The quantum channel N B is said to be more capable than the environment, or just more capable, if the quantum capacity of the complementary channel to the environment is zero, i.e., Q(N E ) = 0. (3)

Definition 2
The quantum channel N B is said to be less noisy than the environment, or just less noisy, if the private capacity of the complementary channel to the environment is zero, i.e., Since C p (N E ) ≥ Q(N E ), less noisy implies more capable. By using the eigenvalue decomposition of ρ A n , we can rewrite the coherent information as where we set |X | = dim H A . Thus, by noting Eq. (2), the quantum channel N B being more capable can be also described as holds for every n ≥ 1. Furthermore, by noting Eq. (1), the quantum channel N B being less noisy can be also described as holds for every n ≥ 1. Eqs. (6) and (7) resemble the definitions of more capable and less noisy for classical channels [19], and it is justified to call quantum channels satisfying Eqs. (3) or (4) more capable or less noisy.
In [19], alternative description of less noisy, less divergence contracting, was introduced, and we can also extend such a description to the quantum channel. The quantum channel N B is said to be less divergence contracting if holds for every n ≥ 1. Later, we will show that the quantum channel is less noisy if and only if less divergence contracting (Proposition 4). This alternative description plays a crucial role when we prove Theorem 2.
The followings are our main results.
Theorem 1 Suppose that the quantum channel N B is more capable. Then, we have Theorem 2 Suppose that the quantum channel N B is less noisy. Then, we have When N B is conjugate degradable, we can show that C p (N E ) = 0 as follows. Suppose that the sender sends a state ρ i that corresponds to the message i, and the eavesdropper 2 uses a POVM {M i }. Then, for the entrywise complex conjugate POVM {M i }, we have , where the last equality follows becauseM T i = M i and (C ⊗n • N ⊗n E (ρ i )) T = N ⊗n E (ρ i ). Thus, the legitimate receiver can get exactly the same information as the eavesdropper and the private information transmission to the eavesdropper is impossible. From this argument, conjugate degradable implies less noisy.
When the quantum capacity of the channel is 0 but it can be used to share bound entanglement, then the channel is called a binding entanlement channel [22]. Especially when the channel produces a positive partial transpose (PPT) bound entanglement, the channel is called PPT binding entanglement channel. If a complementary channel is a binding entanglement channel, then the main channel obviously belongs to the more capable class. Since the complementary channel of the conjugate degradable channel can only produce PPT bipartite state [13], if there exists a conjectured negative partial transpose (NPT) binding entanglement channel, the complementary of such a channel belongs to the more capable class but not to the conjugate degradable class.
It is known that there exists a channel such that the quantum capacity is zero (PPT binding entanglement channel) but the private capacity is strictly positive [18]. Let the complementary channel N E be such a channel. Then the channel N B is more capable but not less noisy 3 . Thus, the more capable class is strictly broader than the less noisy class. However, it is not yet clear whether the less noisy class is strictly broader than the degradable or conjugate degradable classes.

IV. PROOF OF THEOREMS
A. Properties of C In this section, we summarize the properties of C (1) p (N B ) and Q (1) (N B ) when Eqs. (6) or (7) hold for n = 1. For n ≥ 2, we can also show similar properties of C (1) p (N ⊗n B ) and Q (1) (N ⊗n B ) when Eqs. (6) or (7) hold for each n by considering n times extension of N B . The following properties can be regarded as quantum extensions of the properties shown for classical channels in the literatures [2,19,23,24] Proposition 1 Suppose that Eq. (6) holds for n = 1. Then we have be an eigenvalue decomposition. LetX be the random variable on U × X such that where Since P U and {ρ u A } are arbitrary, we have The opposite inequality is obvious from the definitions of C Proposition 2 Suppose that Eq. (6) does not hold for n = 1, and the density operator ρ * A maximizing I c (A B) is full rank. Then, we have Especially when dim H A = 2 and C (1) p (N B ) > 0, the sufficient and required condition on is that Eq. (6) holds for n = 1.
Proof. Since Eq. (6) does not hold for n = 1, there existŝ ρ A such that Since ρ * A is full rank, there exists 0 < λ < 1 such that ρ * A − λρ A is positive semidefinite. We construct P U and {ρ u A } as follows. Let be eigenvalue decompositions. Let U = {0} ∪ X . Then, we set P U (0) = λ, P U (u) = β u for u ∈ X , ρ 0 A =ρ A , and ρ u A = |φ u φ u | for u ∈ X . LetX be the random variable on U × X such that where Eq. (10) follows from I(X; B|U = 0) − I(X; E|U = 0) = I c (A B)ρ < 0 and Eq. (11) follows from Next, we show the latter statement of the proposition. The sufficient condition follows from Proposition 1. Suppose that holds. Since C Then, since dim H A = 2, ρ * A must be full rank. Thus, the required condition follows from the former statement of the proposition.

Proof. Let
be an eigenvalue decomposition.
Then, let U = {1, . . . , m}, P U (u) = p i ,X be the random variable on U × X such that It should be noted that the coherent information is known to be concave if the quantum channel N B is degradable [3, Eq. Proposition 4 The following two conditions are equivalent.
Next, we show that (ii) implies (i). For any P U and {ρ u A }, we have Since Eq. It is a straight forward consequence of Proposition 1. Since N B is more capable, Eq. (6) holds for every n ≥ 1. Thus, we have C (1) p (N ⊗n ) = Q (1) (N ⊗n B ) for every n ≥ 1, and the statement of the theorem follows from Eqs. (1) and (2).

C. Proof of Theorem 2
Since less noisy implies more capable, by Theorem 1, it suffice to show Q(N B ) = Q (1) (N B ). For this purpose, we will show that holds for every n ≥ 1. For any input state ρ A k A ℓ on H ⊗(k+ℓ) A , let ρ A k and ρ A ℓ be the partial traces. Then, we have Since Eq. (7) holds for n = (k + ℓ), by (n times extension of) Proposition 4, Eq. (8) also holds for n = (k+ℓ), which implies Thus, Eq. (14) can be proved by induction.

V. CONCLUSION
In this paper, we introduced two new classes of quantum channels, which we call more capable and less noisy. For the more capable class, we showed that the private capacity and quantum capacity coincide. For the less noisy class, we showed that the private capacity and quantum capacity can be single letter characterized.
Our results shed light on further understanding of the private and quantum capacities of quantum channels. However, the conditions such that a certain channel belongs to the more capable class or the less noisy class are hard to verify in general, and we do not yet know whether there exists a channel that belongs to the less noisy class but not to the degradable or conjugate degradable classes, which is an important future research agenda.


Design and early experience with promoting user-created data in CMS

The Computing Model of the CMS experiment [1] does not address transfering user-created data between different Grid sites. Due to the limited resources of a single site, distribution of individual user-created datasets between sites is crucial to ensure accessibility. In contrast to official datasets, there are no special requirements for user datasets (e.g. concerning data quality). The StoreResults service provides a mechanism to elevate user-created datasets to central bookkeeping ensuring the data quality is the same as an official dataset. This is a prerequisite for further distribution within the CMS dataset infrastructure.

Introduction
In the CMS experiment official datasets and user-created datasets are differentiated. Whereas official datasets are centrally produced by the CMS Data Operations group, users are allowed to produce and store their own datasets containing any kind of data at a Tier 2 Center. In contrast to official datasets, there are no requirements concerning data quality, usefulness or appropriate size to be transfered or stored on tape. User-created data is generally located in the private user space at the their designated Tier 2 and can be registered in a local scope data bookkeeping server, the CMS Dataset Bookkeeping System (DBS) [2]. The provided Grid tools can perform a distributed analysis on user-created data. In principle, this dataset can be analysed by any user of the collaboration, however only at the Tier 2 center hosting the dataset, which naturally has a limited number of job slots. Later the dataset created by the user may become important for many other users or even a whole Analysis Group. To provide better availability it is reasonable to distribute the dataset to additional Tier 2 centers or even to a Tier 1 center for custodial storage on tape. However, the CMS data transfer system (PhEDEx) [3] can only handle official data registered in the central bookkeeping service. Therefore, it is necessary that the user-created dataset becomes an official dataset fitting all the requirements of CMS. The StoreResults service described in this paper provides a mechanism to elevate usercreated datasets ensuring all CMS data quality requirements for official datasets.

The Design of the StoreResults Service
The current system is ad-hoc based around a Savannah request and problem tracker for approvals and on the CMS ProdAgent production framework for distributed processing [4]. The Savannah interface was chosen to be the interim interface to StoreResults, until the new RequestManager service has been put in place. The RequestManager is designed to handle the data processing requests of CMS, such as StoreResults requests. The StoreResults service is implemented as a component of the CMS ProdAgent framework. The structure of its implementation is depicted in figure 2. This ProdAgent component consists of three parts. The RequestQuery part handles the communication with the Request Interface. The requests and their status are periodically synced with the StoreResults table in the ProdAgentDB. New requests are assigned to a physics group Savannah squad according to the given information in the request. These squads are used to notify the physics group conveners or their representatives about new requests awaiting their approval or rejection. Once the request has been approved, a json configuration containing the necessary information is created and the operator is notified that the task is now ready for submission. The actual task submission is manually triggered by the operator and the task is passed to the job handling components of ProdAgent using the ResultsFeeder interface. The bookkeeping, central steering and the monitoring of the processed tasks is done by the StoreResultsAccountant. For this purpose the StoreResultsAccountant periodically queries the internal ProdAgentDB to fetch the current status of the merge jobs processed. In addition, the StoreResultsAccountant also triggers the injection into global DBS and PhEDEx once all merge jobs are successfully done. Afterwards the StoreResultsAccount closes the open request in Savannah by calling a member function of RequestQuery and a notification mail is sent to the requestor, group representatives and the operators.

The Process of Promoting User-created Data
The sequence of promoting user-created data is outlined in figure 3. A user belonging to a specific physics group creates a private dataset or skim using the CMS software tool for distributed analysis (CRAB) [5]. The output of the users jobs is either staged-out to /store/group/<group name> or to /store/user/<user name> at a group or user Tier 2.
Afterwards the dataset is registered in the local scope bookkeeping system, so that it is available for distributed analysis using CRAB.
Since the dataset resides only at one Tier 2 site, the number of analysis slots is limited. To make the dataset available for a larger group of users it is essential to distribute it over several Tier 2 sites. The CMS data transfer systems supports only official datasets registered in global bookkeeping. To proceed the elevation to global bookkeeping the user or a group representative can make a StoreResults request via the Savannah interface. To ensure the usefulness of the dataset an approval by the physics group convener or one of its representatives is required. Once the request has been approved the elevation is started.

Dataset is injected in
Global DBS

Dataset is injected in PhEDEx
Transfer dataset on request by the group Figure 3. The sequence of promoting user-created data in CMS using the StoreResults service.
For data transfer and possibly custodial storage on tape at a Tier 1 center it is mandatory that the individual files of the migrated dataset have a reasonable size, which is recommended to be at least 2 GB. Therefore, merge jobs are submitted to the Tier 2 center, where the user dataset is located. In addition, merging using the CMS software guarantees that the used dataformat is compatible with the Event Data Model (EDM) of CMS. The output of these merge jobs are staged-out to /store/results/<group name> and the dataset becomes the ownership of the requesting group. Since the dataset fulfils now the given requirements concerning file size, usefulness and data quality, it becomes an official CMS dataset and is injected into global bookkeeping and the CMS transfer system PhEDEx. The requesting group is now responsible for initiating data transfers to other locations using PhEDEx and also for the deletion once the dataset becomes obsolete.
The last step of the promotion is the cleanup of the /store/group/<group name> or /store/user/<user name> directory and the invalidation of the dataset in local scope bookkeeping, which is the responsibility of the requesting physics group or user.

Monitoring of the StoreResults Service
To monitor the progress of the elevation a Web Service using CherryPy has been developed. The Web Service is running on each instance of the StoreResults service. The access has been restricted to group representatives for performance and security reasons. The authentication is done by using X.509 grid certificates. A screenshot of the StoreResults monitoring website is shown in figure 4. On one hand the Web Service allows the physics groups to monitor the progress of their requests in real time. On the other hand, it also very useful for the operators of the StoreResults service, to track the handled requests. Since the Standard Output and Standard Error from retrieved jobs are also available through the Web Service, debugging of failed jobs is much simpler. From the technical point of view, the StoreResults monitoring obtains its information about the processed and queued requests by querying the local ProdAgentDB and the bookkeeping database DBS.

Operation of the StoreResults Service
To handle the StoreResults requests, currently two StoreResults instances have been deployed. One is located at RWTH Aachen and the other one at Imperial College London. Since the task submission is done manually, operator actions are required. The StoreResults service is currently operated by three part-time operators. Two people at RWTH Aachen and one at Imperial College London. Furthermore, there are two part-time developers, one at Fermi National Accelerator Laboratory (FNAL) and one at RWTH Aachen.
The primary concern of the physics groups is: "How long does the elevation take?" This question cannot be answered in general, since the time to complete the process depends on various factors. The parent dataset migration from local scope bookkeeping to global bookkeeping takes usually 10 to 60 minutes. This step is necessary to keep the parentage information in the database. A usual task comprises 1-1000 merge jobs. Each jobs lasts between two hours and two days. The elevation task is finished only once the last job is done. However, the merge jobs of StoreResults have a higher priority, since their are running with the CMS Production-Role. Due to recurring Grid errors, each merge job will be resubmitted up to ten times, before it is declared as failed and additional operator action is required. Once all jobs are done, the PhEDEx subscription to a site has to be approved by the operator or the data manager of the Grid site, even if the data resides already at the given site. On average it takes 1-2 days to complete all the steps in elevating a dataset. On average the StoreResults service deals with 1-2 requests per day. Usually, those requests emerge in clusters of 2-10. Both instances handled 333 requests in 2010 so far. About 98% of the request were successfully elevated and the average time is around 46 hours per request.

Future Development
The CMS experiment is going to change its production framework for distributed processing on the Grid to WMAgent, which is an evolution of the currently used ProdAgent system. That means the StoreResults service is going to be completely reimplemented in the WMAgent framework and will also natively support growing datasets.
Accompanying to the transition to WMAgent, the CMS experiment will use the RequestManager to handle data processing requests. The RequestManager combines the possibility of making requests, approvals by authorized people as well as a tracking of their status. Therefore, the currently used Savannah interface will be replaced in the future by an interface integrated in the RequestManager to allow a more convenient handling of StoreResults requests by a unified interface for all data processing requests.


extfieldsofstudy:
Computer Science
Medicine
Engineering

s2fieldsofstudy:
Computer Science
Engineering

text:
Direction-of-Arrival Estimation Based on Joint Sparsity

We present a DOA estimation algorithm, called Joint-Sparse DOA to address the problem of Direction-of-Arrival (DOA) estimation using sensor arrays. Firstly, DOA estimation is cast as the joint-sparse recovery problem. Then, norm is approximated by an arctan function to represent joint sparsity and DOA estimation can be obtained by minimizing the approximate norm. Finally, the minimization problem is solved by a quasi-Newton method to estimate DOA. Simulation results show that our algorithm has some advantages over most existing methods: it needs a small number of snapshots to estimate DOA, while the number of sources need not be known a priori. Besides, it improves the resolution, and it can also handle the coherent sources well.

Introduction
Direction-of-Arrival (DOA) estimation using sensor arrays has been an active research area, playing a fundamental role in many applications involving electromagnetic, acoustic, and communication systems [1]. Many classical algorithms are available, and the popular methods include beamforming [2], MUSIC [3], ESPRIT [4] and the maximum likelihood method [5], etc. The beamforming method has low angle resolution and suffers from the Rayleigh resolution limit. MUSIC, ESPRIT and the maximum likelihood method all rely on the statistical properties of the data, and thus, require a OPEN ACCESS sufficiently large number of samples for accurate estimation. Besides, MUSIC and ESPRIT cannot handle strongly coherent sources, while the maximum likelihood method has high computation costs.
The problem of sparse recovery has evolved rapidly recently [6,7] and it has been applied in DOA estimation with array processing. Gorodnitsky et al. [8] used a weighted least-squares algorithm named FOCUSS for DOA estimation, but this algorithm can only be used for single snapshots. Cotter [9] combined multiple measurement vectors (MMV) and matching pursuit (MP) to solve the joint-sparse recovery problem in DOA estimation, but it has low angle resolution. JLZA-DOA is proposed in [10]; it minimizes a mixed L 2,0 norm to deal with the joint-sparse recovery problem, and a fixed point method is used for DOA estimation. This algorithm doesn't satisfy numerical stability, as matrix inversion is inevitable in every iteration. Stoica et al. [11] presented a novel SParse Iterative Covariance-based Estimation approach, abbreviated SPICE. However, this algorithm needs more snapshots to estimate DOA. Wide-band covariance matrix sparse representation (W-CMSR) is proposed in [12] for DOA estimation of wideband signals. So far, the most successful joint-sparse recovery algorithm for DOA estimation is L1-SVD [13,14]. It combines the SVD step of the subspace algorithms with a sparse recovery method based on l 2,1 -norm minimization. However, the number of sources needs be known a priori.
In this paper, we present Joint-Sparse DOA estimation, abbreviated as JSDOA, for sensor array DOA estimation. First, DOA estimation is cast as a joint-sparse recovery problem. Then, L 2,0 norm is approximated by the arctan function to represent spatial sparsity and DOA estimation can be obtained by minimizing the approximate L 2,0 norm. Finally, the minimization problem is solved by a quasi-Newton method to estimate DOA. The proposed algorithm has some advantages over most existing methods: it needs a small number of snapshots to estimate DOA, an the number of sources need not be known a priori. Besides, it improves the probability of resolution, and it can also handle coherent sources well.
The outline of the paper is as follow. In Section 2, the DOA estimation problem is formulated. The new algorithm, called JSDOA, is proposed in Section 3. In section 4, the validity of the proposed algorithm is proved by a number of simulations. Finally, conclusions are presented in Section 5.

Joint-Sparse Recovery for DOA Estimation Problem
Because sources are sparse in space, DOA estimation with sensor arrays can be expressed as a joint-sparse recovery problem. Let Ω denote the set of possible locations, { } 1 N n n θ = denotes a grid that covers Ω. We assume that the grid is fine enough such that the true location parameters of the existing sources lie on the grid. Let: x t x t x t = Then the received signal y m (t) at the m th sensor can be written as: can be expressed as: when the number of snapshots is denoted as T, Equation (4) can be written as: As we know, X is row-sparse and only K rows have nonzero elements, so DOA estimation can be obtained by a joint-sparse recovery problem, which is also called a multiple measurement vectors (MMV) problem. Using l p,q norm to express joint sparsity, the MMV problem can be converted to l p,q norm minimization: where p and q are non-negative, and , p q X is defined as: Concretely, Eldar and Mishali [15] use p = 2, q = 1. Chen and Huo [16] study for p = 1, q = 1. The above algorithms do not perform very well for ‖X‖ 2,1 and ‖X‖ 1,1 can't sufficiently reflect joint sparsity. Considering ‖X‖ 2,0 can reflect joint sparity sufficiently, we minimize ‖X‖ 2,0 norm to solve the MMV problem. However, ‖X‖ 2,0 norm minimization can hardly be solved directly. Therefore, in this paper we approximate ‖X‖ 2,0 norm by an arctan function and estimate DOA by solving an approximate ‖X‖ 2,0 norm minimization problem.

Joint-Sparse DOA Estimation Algorithm
In this section, the DOA estimation problem is converted into an approximate ‖X‖ 2,0 norm minimization problem. Firstly, the L 2,0 norm is approximated by an arctan function to construct an approximate ‖X‖ 2,0 norm minimization problem. Then this problem is solved by a quasi-Newton method to estimate DOA.

Basic Idea of the Proposed Method
, we have: Considering the following arctan function: where δ is a positive parameter and s is a variable parameter. Then f δ (s) has the following property: From Equations (8) and (10), then: So DOA can be obtained by solving the approximate 2,0 X norm minimization: where δ is a small positive constant. Because noise is unknown, we synchronously hope to minimize F δ (X) and . Then Equation (12) is converted into a multiple objective optimization: Using the linear weighting method, Equation (13) can be written as: where ‖•‖ F denotes the Frobenious norm and the parameter λ will be discussed in Section 3.2.
When the parameter δ is very small, the objective function F δ (X) is highly unsmoothed and contains a lot of local minimization, so its global minimization is not easy. On the other hand, if the parameter δ is larger, the objective function F δ (X) is smoother and contains less local minimization. In order to obtain global minimization, we select a decreasing sequence for δ, denoted as: x δ − is used as the initial value for δ = δ j , thus we hope the proposed algorithm can escape from getting trapped into a local minimum and reach the global minimization for a small value δ = δ J .
For some fixed value δ = δ j , minimization problem Equation (14) is solved by the quasi-Newton method in this paper. One of the most successful quasi-Newton methods is the BFGS algorithm, which is second order convergent and has good numerical stability. Therefore, we use the BFGS algorithm to solve Equation (14) for some fixed value δ = δ j .
The conjugate gradient for a matrix variable is defined as: where X R and X l denote the real part and the imaginary part respectively, X * denotes the conjugate of X.
Then the conjugate gradient of L δ,λ (X) (for its derivation see Appendix I) can be expressed as: Solving by the BFGS algorithm, the main iterative steps are as follows: (1) Search step length t k , satisfy:
(2) Select parameter δ: In this paper, we set , we hope parameter δ 1 satisfies: In order to save computation cost, we set However, if δ J is too small, ( ) J F x δ will be sensitive to noise, so δ J shouldn't be too small.
(3) For coherent sources: Many of popular methods, such as MUSIC and ESPRIT require the assumption that sources are uncorrelated, because the source covariance matrix remains nonsingular so long as none of these sources are coherent. However, the algorithm proposed in this paper doesn't need to compute the source covariance matrix, and joint sparsity is used to estimate DOA, so the proposed algorithm can handle highly coherent sources as well.
(4) Limitation: The algorithm proposed in this paper suffers from two problems. First, there exist no clear-cut guidelines for the selection of δ J , especially when knowledge of noise is unknown. Second, the proposed algorithm can't perform well when the SNR is relatively low.

Simulation Results
In this section, we present several numerical simulation results to illustrate the performance of the proposed algorithm. First, we compare the spatial spectrum of JSDOA to those of beamforming [2], MUSIC [3], and L1-SVD [13] under various snapshot scenartios. Next, we compare the spatial spectrum with more than two sources. Then, the probability of resolution is compared. Finally, we compare the spatial spectrum for coherent sources. In the following numerical simulations, we

Spatial Spectrum Comparison under Various Snapshots
In this simulation, we compare spatial spectrum of different algorithms under various snapshots. We consider two uncorrelated sources located at 13° and 20° with SNR = 10. For JSDOA, we set  Figures 1(a,b), we set T = 10 and T = 5, respectively.
It can be seen from Figure 1 that JSDOA and L1-SVD can resolve the two sources as T = 10. However, only JSDOA can resolve the two sources as T = 5.

Spatial Spectrum Comparison with More Than Two Sources
In this simulation, we compare the spatial spectrum of different algorithms with more than two uncorrelated sources. Set SNR = 10, T = 20,

Probability of Resolution Comparison under Various Conditions
In this simulation, we compare probability of resolution for two uncorrelated sources. JSDOA, MUSIC and L1-SVD are used for comparison and we set T = 5,  Figure 3(a), one source is fixed at 10°, the other source is located at 10° + Δθ. vary from 0 to 30. In Figure 3(b), two sources are located at 10° and 20°. SNR vary from −5 to 20. It is shown from Figure 3 that JSDOA has the higher probability of resolution than MUSIC and L1-SVD.

Spatial Spectrum Comparison for Coherent Sources
In this simulation, we compare spatial spectrum for completely coherent sources. Set SNR = 10, T = 20,  In Figure 4(b) three coherent sources, located at 10° and 20°, are considered. It is shown from Figure 4 that JSDOA and L1-SVD can be used to estimate DOA for coherent sources and that JSDOA has higher resolution than L1-SVD.

Conclusions
In this paper, a joint-sparse recovery algorithm, called JSDOA, is proposed to estimate DOA with sensor arrays. We pose the DOA estimation problem as a joint-sparse recovery problem, which can be solved by minimizing the approximate L 2,0 norm. In particular, the L 2,0 norm is approximated by an arctan function to represent spatial sparsity and the approximate L 2,0 norm minimization problem is solved by a quasi-Newton method. Finally, the proposed algorithm is examined by simulations. Several advantages over existing DOA estimation algorithms were identified. It can perform well with a limited number of snapshots, while the number of sources need not be known a priori. Besides, it improves the resolution, and it can also handle the coherent sources well.

Acknowledgment
This work was supported by the National Natural Science Foundation (No. 61072120), the program for New Century Excellent Talents in University (NCET).


A Distributed Agent Based System to Control and Coordinate Large Scale Data Transfers

We present a distributed agent based system used to monitor, configure and control complex, large scale data transfers in the Wide Area Network. The Localhost Information Service Agent (LISA) is a lightweight dynamic service that provides complete system and applications monitoring, is capable to dynamically configure system parameters and can help in optimizing distributed applications. As part of the MonALISA (Monitoring Agents in A Large Integrated Services Architecture) system, LISA is an end host agent capable to collect any type of monitoring information, to distribute them, and to take actions based on local or global decision units. The system has been used for the Bandwidth Challenge at Supercomputing 2006 to coordinate global large scale data transfers using Fast Data Transfer (FDT) application between hundreds of servers distributed on major Grid sites involved in processing High Energy Physics data for the future Large Hadron Collider experiments.

INTRODUCTION
The High Energy Physics community engaged in CERN's Large Hadron Collider (LHC) is preparing to conduct a new round of experiments to probe the fundamental nature of matter and space-time, and to understand the composition and early history of the universe. The decade-long construction phase of the accelerator and associated experiments is now approaching completion, and the design and development of the computing facilities and software is well-underway. The LHC is expected to begin operations in 2007. The experiments face unprecedented engineering challenges due to the volume and complexity of the experimental data, and the need for collaboration among scientists located around the world. The massive, globally distributed datasets which will be acquired, processed, distributed and analyzed are expected to grow to the 100 Petabyte level and beyond by 2010. Distribution of these datasets will require network speeds of around 10-100 gigabits per second (Gbps) and above. The data volumes are expected to rise to the Exabyte range, and the corresponding network throughputs to the 100 Gbps -1 Terabit/sec range, by approximately 2015. In response to these challenges, the Grid-based infrastructures developed by collaborations in the US, Europe and Asia such as EGEE, OSG and Grid3 provide massive computing and storage resources. However, efficient use of these resources is tightly connected with the efficient use of the deployed networking equipment.
Even with high-speed networking resources deployed between hosts involved in the LHC experiments, an important problem that needs to be solved is related to End-to-End achievement. During the "High Speed Data Gathering, Distribution and Analysis for Physics Discoveries at the Large Hadron Collider" demonstration at SuperComputing 2006 (SC06) Bandwidth Challenge (BWC) our agent based solution provided the control plane that helped in the achievement of new records for sustained data transfer.

THE ARCHITECTURE OF LISA
LISA (Localhost Information Service Agent) is a lightweight framework that can help in optimizing other applications by means of monitoring services. It consists of a number of service agents that can be easily deployed on any workstation, independent of the local architecture or operating system. The LISA agent automatically detects the architecture on which it is deployed and is capable to dynamically load the binary modules necessary to perform monitoring services. LISA is based on Java technologies, any agent being able to operate independent of the operating system platform. The components of the framework are presented in figure 1.
The service agent is composed of a core system and a set of deployable modules. The core system is responsible for managing monitoring modules in the framework. It dynamically discovers available modules and can load or unload them at any time. The modules can be loaded or unloaded without the need of restarting the system and without affecting the rest of the running modules. Another interesting feature of the engine is that is capable of sensing when parameters are updated. When such updates occur the engine is capable of triggering sets of actions such as reloading the modules or updating the running parameters.
The core system is capable of working with remote repositories of modules. From remote locations it can dynamically trigger updates of module. The core system also provides the capability of integrating the agent with another external application using sets of exported API methods. Every module in the system can export such functionality to any external application.  Another aspect of the implementation of the core system is that it is permanently monitoring itself. This capability makes LISA capable of detecting potential problems and choosing from a number of possible solutions to counteract them. For example if a loaded module stops due to some running conditions or programming fault the core system is capable of dynamically restarting it. The core system is even capable of restarting the agent altogether, for example when the Java virtual machine encounters a problem such as the available memory wears out or the running threads come to a deadlock condition.
The LISA Service Agent incorporate a number of monitoring and controlling modules that can be used to completely supervise the end-user system, the monitoring capabilities ranging from the hardware architecture and the running parameters of the workstation to local networking parameters. New modules can be easily constructed and added by users and this capability provides increased monitoring functionality to the system.
One of the already implemented and available modules offers complete system monitoring functionality. This module can monitor environmental information such as the type of the operating system and version, the user under which LISA runs, JVM version, local IP address or, even better, if the address is a private one it can retrieve  the public address of the access point for the local network, and the AS to which the station belongs. This module is also able to monitor the health of the LISA agent using JMX technology. Figure 2 presents the type of information that this module is capable of gathering from the workstation where LISA is deployed.
One other available monitoring module supervises the status of the workstation health parameters. The monitored information is classified into information regarding the CPU utilization (how much of the CPU is used by the user applications, by the system, how much is idle, etc.), regarding memory utilization (free memory and total available memory in the system, page swapping), regarding disk usage (free disk space and total available disk space), regarding the current load of the system (load1, load5, load15) or the number of running processes and information regarding networking (in and out traffic for each available network interface). This module is dynamically adjusting itself based on the architecture of the system counters of the underlying workstation on which LISA is running. Examples of the capabilities offered by this module are presented in figure 4. LISA also comes with a module that monitors the hardware configuration of the local workstation. The information collected by this module is independent of the underlying operating system offers. Figure 5 presents an example of the kind of ability presented by this module.
One other module is capable of monitoring and adjusts the networking performances of the local station. This module reads information regarding the network interfaces available in the system such as traffic patterns, addresses and capacities and current settings. Also it monitors specific parameters related to the installed TCP/IP stack. It reads parameters regarding the health of the IP protocol (packets sent and received, bad packets, resent packets, etc.), the health of the TCP protocol (sent and received segments, fragmentation, etc.) and the health of the UDP protocol. This module is also able to adjust on user request some of the parameters that are used for tuning the TCP/IP stack in order to optimize the networking performances of the local workstation. Examples of the capabilities offered by this module are presented in the figure 7.  One other implemented module can be used as a shell interface between any user having correctly defined rights and the local workstation on which the LISA agent is deployed. Using this module the user, meaning the user running the GUI client or an external application, can execute commands and read the returning outputs. And another module is able to test the available networking throughput between the local workstation and any other remote station, using IPERF. Figure 6 presents results obtained using this monitoring module.
All these modules are meant to make LISA an application which best describes the local station in order to help other applications optimizing performances by means of monitoring services. There are actually two defined ways in which an application can integrate the LISA services. Any application can use the API functions that each module is exporting. In this way any application can be tightly coupled with the LISA framework. Another modality is that the application connects itself to the service daemon. In fact this is how LISA is connecting to the higher level MonALISA monitoring services and also this is how the GUI client connects to the LISA services. In both approaches the LISA agent is sending to the registered interfaces the monitoring results and can receives commands to trigger internal actions.
The lightweight monitoring framework also consists of a GUI client which can be used by any user to keep track of the current monitoring results. The GUI client is using Java WebStart technology and because of that it can be executed on any local station.
The GUI client can connect to any number of LISA services using networking communication. In this way the status of the services is monitored from any remote location. When connecting to a remote monitoring agent the user running the GUI client have to authenticate in order to provide a set of rights. Based on this access rights the user can have only the capability to receive results or it can more advance rights such as to modify running parameters or rights to start or stop modules.
Each module loaded in the LISA framework is able to define a graphical interface. The client does not start with any of the graphical interfaces defined by any of the running modules. All these interfaces are dynamically loaded into the GUI client. The interfaces can be loaded from any remote locations, from repositories of modules or even can be loaded from the agent itself where a particular module is running. When the user connects to a running LISA agent it will not receive from the beginning all monitored results. When the user becomes interested in a particular module (by clicking on that particular module) the dynamic interface is loaded in the GUI system and the client starts receiving monitoring results.
The client is able to preserve the list of known LISA agents from one run to another. This combined with the aspect that the GUI client is able to remember certain repeated actions means that the system is designed to be user-friendly.
Because LISA is a part of MonALISA there is another module which sends the monitored values back to one or more MonALISA service. For that this modules uses the ApMON module provided in the MonALISA framework. The MonALISA can help in the way that the monitored result values can be processed at a higher level using the MonALISA inner functions.

FAST DATA TRANSFER ARCHITECTURE
FDT is an Application for Efficient Data Transfers which is capable of reading and writing at disk speed over wide area networks (with standard TCP). It is written in Java, runs an all major platforms and it is easy to use. FDT is based on an asynchronous, flexible multithreaded system and is using the capabilities of the Java NIO libraries. FDT performs the streaming of a dataset (list of files) continuously, using a managed pool of buffers through one or more TCP sockets. It uses independent threads to read and write on each physical device and transfers data in parallel on multiple TCP streams, when necessary. It uses appropriate-sized buffers for disk I/O and for the network FDT restores the files from buffers asynchronously and resumes a file transfer session without loss, when needed.
FDT can be used to stream a large set of files across the network, so that a large dataset composed of thousands of files can be sent or received at full speed, without the network transfer restarting between files.
The architecture of the Fast Data Transfer application is described in figure 8.

LISA AT SC06 BANDWIDTH CHALLENGE
The future of Grid computing lies in the ability to collaborate and transfer petabytes of information around the world. Experiments such as CMS and Atlas are currently constructing detectors for CERN's Large Hadron Collider (LHC) which will unite thousands of physicists worldwide. In order to support such a global project, worldwide interconnections between large databases, mainframes servers, clusters and terminals are required to provide access, processing and analysis of the petabyte sized datasets that experiments such as LHC will provide. The interconnections between machines need to be capable of transferring the traffic that the new datasets will generate.
To keep this flow of information free from disruption, it is not only the servers and clusters that have to be kept operational and optimal, but also the underlying networks that feed the processing farms. It is therefore critical to understand the performance of the underlying networks in order to plan, prepare and better utilize the network. The Fast Data Transfer (FDT) application was designed to meet this scope, to optimize long-range network transfers. During the "High Speed Data Gathering, Distribution and Analysis for Physics Discoveries at the Large Hadron Collider" demonstration at Fig. 9. Results obtained by consecutively adding one more data stream to the overall network transfer flux.
SuperComputing 2006 (SC06) Bandwidth Challenge (BWC), an international team of physicists, computer scientists, and network engineers led by the California Institute of Technology, CERN, and the University of Michigan joined forces to set new records for sustained data transfer between storage systems. Using FDT for its networking support and LISA as the control plane of the experiment the team achieved a peak throughput of 17.77 Gbps between clusters of servers located at the show floor and at Caltech.
Following the rules set for the SC06 Bandwidth Challenge, the team used a single 10-Gbps link that carried data in both directions. FDT provided sustained total throughput of around 17 Gbps for disk to disk transfer using 10 pairs of small servers (each having 4 SATA HD configured in software raid0 and 1Gb/s network interface) in both directions. The official results obtained are presented in figure 10.
During this event LISA played an important role in establishing the record. Through continuous monitoring and supervising tasks LISA provided the control plane for the challenge. This is a perfect example of how LISA can help optimize other applications.
The networking module was used to continuously  For this event all monitoring module provided precious information about the state of the system. The parameters reported by the host monitoring module provided good insights of the internal state of the computing nodes between which the transfers took placed. The bash module was used to execute diverse commands on the remote nodes in a simple and intuitive way.
LISA also provided one module that was used in close relation with the FDT application. This module reported the state of the application and information pertaining to the network transfers in progress. LISA provided feedback in real-time of the status of the controlled system. LISA was used to dynamically adjust the running parameters of the FDT application, was used to restart the running networking transfers or to change the files being transferred. And LISA performed everything in real time. Operations such as kernel parameters retrieval or asymmetric MTU detection performed by LISA also proved to be crucial in establishing the new records.
The high-energy physics team also carried out several other demonstrations, making good use of the ten wide-area network links connected to the Caltech/CERN booth. Some results are represented in figures 9 and 11. The challenge demonstrated a realistic, worldwide deployment for distributed, data-intensive applications capable to effectively use and coordinate the network resources. In the overall efforts for obtaining the results an important role was played by LISA.

LISA AS PART OF THE EVO SYSTEM
Another example of the way LISA can help other applications in terms of optimization is the relationship between the framework and the EVO videoconference system. As mentioned in (Dobre, et al., 2004), LISA was made part of the functioning EVO system, assisting the EVO peers to dynamically detect the best reflectors to which to connect to. The best reflectors are read from MonALISA repositories and are updated in time so that the choosing is performed only from the available reflectors. Then the best reflectors are chosen based on their network location (Network domain, AS domain, Country, Continent) and on their current load values, number of currently connected clients, current network traffic. This insures a load balancing order in choosing the reflectors. Based on this from all the available reflectors the module that chooses some and performs RTT values measurements. In the end the module is able to inform some application (e.g. VRVS client) to reconnect to the most appropriate module if network conditions change and the quality is affected.

CONCLUSION AND FUTURE PLANS
The distributed agent system based on LISA provides the functionality to monitor and dynamically control large scale systems. It was been successfully used to coordinate and optimize large scale data transfers. It is currently used in the EVO videoconference system to optimize the connectivity and performance from the user's perspective. LISA can be used together with the MonALISA framework to provide optimized workflows in distributed grid systems. We are working to provide support for dynamical optical path allocation between end systems.


s2fieldsofstudy:
Computer Science
Business

text:
Cooperative and Distributed Localization for Wireless Sensor Networks in Multipath Environments

We consider the problem of sensor localization in a wireless network in a multipath environment, where time and angle of arrival information are available at each sensor. We propose a distributed algorithm based on belief propagation, which allows sensors to cooperatively self-localize with respect to one single anchor in a multihop network. The algorithm has low overhead and is scalable. Simulations show that although the network is loopy, the proposed algorithm converges, and achieves good localization accuracy.

I. INTRODUCTION
A wireless sensor network (WSN) consists of many devices (or nodes) capable of onboard sensing, computing and communications. WSNs are used in industrial and commercial applications, such as environmental monitoring and pollution detection, control of industrial machines and home appliances, event detection, and object tracking [1]- [3]. In most applications, the data collected by the sensor nodes can only be meaningfully interpreted if it is correlated with the location of the corresponding sensors. The Global Positioning System (GPS) is widely used for localization in outdoor environments [4]. However, GPS is a costly option and is not suitable for power-limited sensor nodes in WSNs. Furthermore, GPS signals do not penetrate well to indoor environments. Therefore, alternatives to GPS localization have been widely studied [4], [5]. In a WSN, nodes whose positions are known are called "anchors". By making use of pairwise range or angle measurements between anchors and/or other sensor nodes whose positions are unknown, sensors with no access to GPS can perform self-localization. Typical techniques include the use of time-of-arrival (TOA), time-difference-of-arrival (TDOA), received-signal-strength (RSS) and/or angle-of-arrival (AOA) information in triangulating the location of a node. This is usually studied in line-of-sight (LOS) environments [6], [7].
However, LOS signals do not always exist in urban or cluttered environments, where signals usually experience multiple refections and diffractions. Such signals are referred to as nonline-of-sight (NLOS) signals and are commonly encountered in both indoor (e.g., residential buildings, offices and shopping malls) and outdoor (e.g., metropolitan and urban) environments. NLOS  extensively investigated [8]- [14], but most of the algorithms in the literature focus on locating one single sensor with several anchors. Since they require each sensor to have direct signal paths to anchors in the network, such algorithms cannot be applied in network-wide localization. On the other hand, several algorithms for network-wide localization have been proposed in the literature [5], [15]- [18]. Unfortunately, only LOS signals are considered in these algorithms.
Distributed localization algorithms for multipath environments were proposed in [19], [20], where NLOS error is modeled as a positive bias in range and angle measurements, and its statistical characteristics are inferred by numerical methods, such as bootstrap sampling in [19], and particle filters in [20]. One of the major disadvantages is that these Bayesian inference techniques require a large number of observations and are computationally expensive. Generally, the statistical model for NLOS errors depends on various factors, such as signal bandwidth, propagation medium and environment temperature. Different models including the uniform, exponential and Rayleigh distributions, have been proposed in the literature [21], [22].
Instead of modeling NLOS errors in multipath environments as random biases, geometric analysis can be applied in pairwise localization [12]- [14], where measurements of different paths are modeled using closed-form expressions, which significantly simplifies the system model. In this paper, we derive a distributed localization algorithm based on range and direction measurements at each node, and where nodes exchange information to cooperatively perform self-localization relative to a fixed reference node. We show through simulation that by exchanging limited information, all the nodes in the network can perform localization to a good accuracy. We compare the performance to that achieved without cooperation, and show that cooperation amongst neighboring nodes significantly improves the localization accuracy.
The rest of this paper is organized as follows. In Section II, we define the system model. We describe our algorithm in Section III, and provide simulation results in Section IV. In Section V, we summarize and conclude. The objective of each S i is to perform self-localization relative to S 0 .
In the following, we consider two nodes S i and S j . Similar to [14], we describe a model to relate the range and direction measurements at each node to their positions. Suppose that there are R LOS or NLOS paths between S i and S j . An example of a single-bounce scattering path is shown in Figure  1, where the signal from S j to S i is reflected at a nearby scatter, and the communication link between two nodes is assumed to be symmetric. Let d r ji be the distance measured by S i using the time-of-arrival information of the signal along the r th path from S j , and θ r ji be the corresponding angle-ofarrival information. Nodes S i and S j exchange measurements with each other, so that both nodes have the measurements {d r ij , d r ji , θ r ij , θ r ji } R r=1 . Consider the r th path between S i and S j . Given the position of S j and {d r ij , d r ji , θ r ij , θ r ji }, the position of S i cannot be determined with certainty even in the absence of measurement and communication noise. As shown in Figure 1, the estimated position for S i can be any point along the line AB. If there are multiple paths between S i and S j from non-parallel scatters, the position of S i can be found as the intersection point of two such lines. Suppose that there is no measurement noise, then a straightforward geometric consideration shows that where p A and p B are the positions of A and B respectively. A vector perpendicular to p A − p B , is Since both A and S i are on the line AB, we have n T , which can be further simplified as where We have made use of the fact that d r ij = d r ji for symmetric communication links between S i and S j in (3). We have not factored in measurement and communication noise up to this point. Let the corresponding noisy measurements be {d r ji ,θ r ji ,θ r ij }. Modeling the total effect of noise as a Gaussian random error r ji , we havẽ We assume that measurements for the R paths are such that g θr ji ,θ r ij : r = 1, . . . , R are linearly independently, otherwise some of the paths are duplicates of each other. We also assume that the noise terms r ji are i.i.d. Gaussian random variables with zero mean and variance σ 2 . Stacking the measurements from all R paths into a vector where G † ji is the Moore-Penrose pseudo-inverse of the matrix G ji . If there are more than one signal paths between nodes S i and S j , we have The posterior distribution of the node locations is given by [14] proposes a least square estimator to localize a single node with respect to a reference node. However, to localize every node in a network, it is necessary to consider the interaction between S i and all the other nodes {S j } M j=0,j =i . The MAP estimator for s i is given by (7) on top of next page. The joint posterior distribution in (7) depends on interactions amongst all the variables and is difficult to compute by brute-force integration.
In order to provide a computationally efficient algorithm to calculate the marginal distributions, we observe that the state at any sensor depends directly only on its neighboring sensors whose number is usually far less than the total number of variables. To explore such conditional independence structure, we make use of belief propagation in the following and propose an efficient algorithm which computes a set of marginal distributions from the joint posterior distribution without performing a full integration as in (7).
Remark 1: When there exists a LOS path between S i and S j , it is easy to see that |θ r ji − θ r ij | = π and d r ji = |s i − s j |, which is a special case of (3).
Remark 2: When there exists paths with multiple bounces between S i and S j , a two-step proximity detection scheme suggested in [14] can be applied to detect and discard such paths, leaving measurements from either LOS and/or singlebounce paths for processing.

PROPAGATION
We use Belief Propagation (BP) on a factor graph [23] in this paper. An example for a network with 5 nodes is shown in Figure 2. Each random variable s i is represented by a variable node (circle). The interaction between two sensors S i and S j is represented by a factor node (square) connected to both variable nodes s i and s j . We split the interaction between S i and S j into two factors f ji Since messages only flow in one direction along the edges in the factor graph, they can be broadcast by the sensor nodes.
Without loss of generality, suppose the sensor S 0 is the anchor with a known position (0, 0). For each variable s i , i ∈ {1, · · · , M }, the marginal posterior distribution m(s i ) is found by iterative belief propagation, where two kinds of messages are involved, where B i is the index set of S i 's neighboring sensors.
Therefore, setting the initial belief b i (s i ) to be the corresponding prior distribution p(s i ), beliefs (8) and messages (9) are iteratively updated at each sensor, and the estimation for s i is found by maximizing the converged belief b

A. Derivation for closed-form beliefs and messages
To derive closed-form expressions for (8) and (9), we first consider the message from the anchor S 0 to a neighboring sensor S i . Since s 0 ≡ (0, 0), its belief is a constant and can be represented as b (l) 0 (s 0 ) = δ{s 0 , (0, 0)}, where δ{·, ·} is the Krocnecker delta function. Therefore, the message from S 0 to S i will be constant over all iterations, and is given by where ν 0i = G † 0i d 0i and W 0i = σ 2 Σ 0i with Σ 0i = G † 0i (G † 0i ) T . Other the other hand, for all nodes S j other than S 0 , we set the initial belief of S j as a Gaussian distribution with zero mean and large variance, i.e., for j = 1, · · · , M , b where α is a large positive value of at least an order of magnitude larger than the dimension of the environment in which the sensor nodes are located. Let λ max ij be the largest eigenvalue of Σ ij . We assume that α ≥ 2σ 2 max i,j λ max ij . As the BP algorithm proceeds, the belief b Algorithm 1 Distributed Localization in Multi-path Environments 1: Initialization: 2: Set the position at the anchor S 0 as s 0 = (0, 0). 3: Set µ (0) i = (0, 0) and where ν  The belief of s i after the l th iteration then follows from (8). Since all the messages in the product in (8) are Gaussian distributions, we have b and µ (l) At the end of the l th iteration, sensor S i estimates its position by maximizing the belief b i (s i ) is a Gaussian distribution, the estimator for s i at the l th iteration is given byŝ i . This iterative procedure is formally given in Algorithm 1. Notice that factor nodes are virtual, and are introduced only to facilitate the derivation of the algorithm. In practice, the update is performed at each individual sensor directly.

IV. SIMULATION RESULTS
Numerical simulations are conducted to validate the effectiveness of our proposed algorithm. We consider a network   First, we consider scenarios where scatters are orthogonal. We compare the performance of cooperative and pairwise localization. The cumulative distribution functions for absolute localization error of each sensor are shown in Figure 3. Corresponding to the factor graph in Figure 2, sensors S 1 and S 2 are directly connected to the anchor. Sensors S 3 and S 4 do not have any paths to the anchor. Instead each has a NLOS path to S 1 or S 2 respectively, and a NLOS path between themselves. In pairwise localization, S 3 localizes using only measurements from S 1 , and S 4 localizes with respect to S 2 . It can be seen from Figure 3(a) that S 3 and S 4 are localized with larger errors than S 1 and S 2 , and this is because errors are accumulated over hops. In cooperative localization, S 3 and S 4 exchange information and incorporate measurements from the NLOS path between themselves. As shown in Figure 3(a), the cooperative localization achieves better performances with more than 90% of the localization errors less than 2m and all errors smaller than 3m. Moreover, results for estimation on y-coordinates are shown in Figure 3(b). It can be seen that both schemes have similar performances for sensors directly connected to the anchor (e.g., S 1 and S 2 ), and cooperation localization for S 3 and S 4 achieves better performances. The mean absolute errors for both schemes are further shown in Table I. Second, we consider scenarios where scatters are horizontal or at angle 45 • to the horizontal. As can be seen from Figure 4, cooperation among neighboring sensors improves performance on both xand ycoordinates. We also note that compared with Figure 3(b), estimation errors for y-coordinates deteriorate due to correlation between measurements on the vertical direction. On the other hand, we investigate the convergence rate for our proposed algorithm in this biorthogonal scenarios, and the result is shown in Figure 5. Simulations are also conducted when scatters are at 10 • , 20 • and 30 • to the horizontal. Similar results as in Figure 5 are obtained and hence omitted here. These numerical results suggest that even for scenarios with non-orthogonal scatters, the mean of belief at each sensor converges in the proposed algorithm.

V. CONCLUSIONS
In this paper, we propose a distributed algorithm based on belief propagation for network-wide localization in multipath environments. The proposed algorithm requires communications only between neighboring sensors, with each sensor processing only information local to itself. The proposed algorithm has low overhead and can achieve robust and scalable localization. By utilizing both TOA and AOA information of the single-bounce scattering paths, we require only one anchor in the whole network, and sensors that do not have LOS/NLOS paths to the anchor can be localized by cooperation with its neighboring sensors. Simulation results show that our proposed algorithm achieve accuracy less than 1m in a 10m×10m square area.


extfieldsofstudy:
Physics
Computer Science

text:
Parallelization of maximum likelihood fits with OpenMP and CUDA

Data analyses based on maximum likelihood fits are commonly used in the high energy physics community for fitting statistical models to data samples. This technique requires the numerical minimization of the negative log-likelihood function. MINUIT is the most common package used for this purpose in the high energy physics community. The main algorithm in this package, MIGRAD, searches the minimum by using the gradient information. The procedure requires several evaluations of the function, depending on the number of free parameters and their initial values. The whole procedure can be very CPU-time consuming in case of complex functions, with several free parameters, many independent variables and large data samples. Therefore, it becomes particularly important to speed-up the evaluation of the negative log-likelihood function. In this paper we present an algorithm and its implementation which benefits from data vectorization and parallelization (based on OpenMP) and which was also ported to Graphics Processing Units using CUDA.

Introduction
The maximum likelihood (ML) fitting procedure is a popular statistical technique used to estimate parameters of a statistical model on a given data sample [1]. Data samples are a collection of N independent events, an event being the measurement of a set of variableŝ x = (x 1 , . . . , x n ) (energies, masses, spatial and angular variables...) recorded in a brief span of time by physics detectors. The events can be classified in different species, which are generally denoted with signals, for the events of interest for their physics phenomena, and backgrounds, all that remains. Each variable x j is distributed for the given species s with a probability distribution function (PDF) P j s (x j ;θ j s ), whereθ j s are free (not constant) parameters of the PDF. If the variables are uncorrelated each other, then the total PDF for the species s is expressed by P s (x;θ s ) = j P j s (x j ;θ j s ).
The ML technique allows to estimate the values of the free parameters, as well the number of events belonging to each species n s , by maximizing the function which is called the extended likelihood function. We should underline thatx i are measured and the P j s functions are well-known, so L only depends on the free parameters we want to fit on the data sample.
The search for the maximum for L can be carried out numerically. Usually, it is used to minimize the equivalent function − ln(L), the negative log-likelihood (NLL). So the NLL to be minimized has the form 1 : that is a sum of logarithms. The most common method used in the high energy physics (HEP) community for the minimization is based on the MIGRAD algorithm inside the MINUIT package. MIGRAD performs the minimization using the variable metric method [2]. This method involves the calculation of the derivatives of the NLL for each free parameter. Since very often we deal with minimizing a function for which no derivatives are provided, MIGRAD is able to estimate the derivatives of the function by finite differences [3]. The whole procedure of minimization requires several evaluations of the NLL, which requires themselves the calculation of the corresponding PDFs for each variable and each event of the data sample. Hence, depending on the complexity of the NLL function with several free parameters, many independent variables and large data samples, the minimization procedure can be very time-consuming. In this case it is important (or even mandatory) to speed-up the evaluation of the NLL [4]. The common software used in HEP community for the evaluation of the NLL is RooFit [5], which is part of the general data analysis framework ROOT [6]. Currently RooFit implements an algorithm for the NLL evaluation which cannot take full advantage of data vectorization and other code optimizations (like function inlining) due to its implementation based on C++ virtual methods [4]. To overcome these limitations, we have designed and implemented a new optimized algorithm, and parallelized it using a data parallelism paradigm implemented with OpenMP. The algorithm has been also implemented to run on a Graphics Processing Unit (GPU) device by using the CUDA language provided by NVIDIA. In this paper we describe the algorithm and the two implementations. Finally we show the comparison of the performance between our CPU implementation and RooFit and the comparison of the performance between CPU and GPU implementations.

NLL Evaluation
The RooFit package is formed of a set of C++ classes constructed on top of the ROOT framework dedicated to likelihood-based analyses. Basically for each mathematical concept there is a corresponding C++ class, e.g. classes for the PDFs and variables definition. Then there is a special class which takes care of finalizing the NLL calculation. Furthermore RooFit provides an interface to the MINUIT package. We should underline that all floating point operations are performed in double precision. Data is organized in memory like a matrix where the columns contain the values for each variable, and the rows represent the values of the variables belonging to each event. All classes for PDFs inherit from a common abstract class, which provides the common interface. So each class has a virtual method to get the value of the PDF. Combinations of PDFs are possible with classes for adding, multiplying and convoluting basic PDFs.
In order to calculate the NLL from the formula (3), the current available RooFit algorithm consists of the following steps (in order): • For a given set of values of NLL free parameters, loop over the events i = 1...N : -read the values of the variables for event i; 1 We omit the N ! term in the expression, which does not depend on the parameters. -calculate the PDFs for the event i; -combine, by means of addition and multiplication, the results of the individual PDFs to calculate the total PDF value for the event i; -calculate the logarithm of the total PDF value, which is the term in the sum of the NLL; -accumulate the terms of the sum.
• Finalize the calculation of the NLL.
The key part of this procedure is the calculation of all PDFs for each event, and then there is a single loop over all events. Since this is done by having recourse to calls of the virtual method of each PDF, this algorithm does not allow particular code optimization, like inlining and data vectorization, and it introduces the obvious overhead due to the virtual method calls.
In order to advantage from code optimization, we redesigned the algorithm to reduce the number of calls to virtual methods. Furthermore, the data is stored differently: the values of each variable are organized in independent arrays, so that we can profit from the coalescing of memory accesses for each variable. The new algorithm follows a different procedure with respect to the RooFit algorithm described above: • For a given set of values of the parameters and a given PDF, we evaluate the PDF on each event of the data sample (which means calculating the PDF on the corresponding arrays of variables), and we save the results of this calculation in an array. So we do a loop over all events i = 1...N and calculate the PDF for each of them. • Repeat the previous step for all PDFs, so we end up with several arrays of partial results (an array for each PDF). Each array of results is composed by N elements, i.e. a result for each event. • Combine, by means of addition and multiplication, all arrays of partial results, corresponding to each event, providing a final array of results, i.e. the array of results of the total PDF. • Calculate the logarithm of the total PDF results.
• Do the sum of the total PDF results and finalize the calculation of the NLL.
The key part of this procedure is the calculation of each PDF for all events, so that instead of one single global loop over the events, now we have independent local loops for each PDF (and their combinations). For the implementation of this new algorithm in RooFit, we add a new virtual method for each PDF class with a reference to the data sample as parameter. Inside this method we perform the local loop over all values of the variables of the corresponding PDF, storing the results of the calculations in an array of partial results. Then the method returns a reference to this array. Since this new virtual method is called just once per each PDF during an NLL evaluation, and then within local loops we perform the calculations of the mathematical functions for all events, we can conclude that the number of calls to virtual methods does not depend by the number of events. Furthermore, thanks to the new data structure organized as arrays for each variable, this code can easily be vectorized for the calculation of each PDF. The loop over the final results of the total PDF to calculate and finalize the NLL evaluation is done in the usual class for the NLL finalization. We should note that a drawback of this new algorithm is that we have to manage all the arrays for the temporary results.

Parallelization Strategy and Implementation
We parallelize the implementation of the new algorithm using OpenMP. Since the iterations in the loops inside the new virtual methods for the calculation of the PDFs are independent, we can parallelize them for the execution on CPUs via the #pragma omp parallel for directive. Each iteration in these loops calculates the corresponding PDFs for a single event. The scheduling of the iterations is static, and each thread does a fixed number of iterations and accesses to consecutive elements of the arrays of input variables and partial results, allowing data vectorization of the loops. Arrays of input variables and results are shared among the threads, so that there is a negligible increment in the global memory footprint of the application when running it in parallel. We parallelize also the sum over the final results of the total PDF inside the class for the NLL finalization (parallel reduction). To mitigate rounding problems in the parallel reduction, which depend on the number of threads, we adopt the algorithm described in [7]. Together with the implementation based on OpenMP parallelization, we also implement a parallelization for a GPU device based on the CUDA language. Basically, the loops inside the new virtual methods for the calculation of the PDFs are substituted by calls to CUDA kernels. Hence, the kernels execute the calculation of the PDFs in parallel tasks, where a task represents the calculation of a given PDF on a single event. Each task gets evaluated by a CUDA thread, i.e. there is a one to one correspondence between tasks and threads. Of course the occupancy of the resources used by each CUDA kernel depends on the corresponding PDF function to be evaluated, which means that an optimization of the occupancy would require a custom task partitioning for each kernel. However, we have checked in simple cases (a kernel used for a Gaussian PDF) that the benefit from such an optimization is small, while it would require more efforts for implementing it. So we prefer to use the current configuration, leaving to next versions of the code the possibility to consider a better partitioning of the tasks per CUDA kernel. Threads are synchronized after the corresponding kernel call for each PDF. All data and calculations execute on the GPU are, as before, in double precision floating point. Arrays of the input variables are copied from the host to the device global memory using synchronous functions. These arrays are read-only during the entire execution of the application, so we can copy them once to the device memory at the beginning of the minimization procedure and then use them for all NLL evaluations. The arrays of partial results for each PDF are only allocated and kept resident in the global memory of the device, i.e. no copy from device to host, except for the array of the final results which has to be copied to the host memory for the final sum of the NLL. Also in this case the communication is done by using synchronous functions. Then, we perform the parallel reduction of the final results using the OpenMP implementation. To summarize, during the minimization there is a single copy of the arrays of the input variables and a copy of the the array of the final results for each NLL evaluation. With this implementation we are able to strongly reduce the time spent for the communication between host and device memories. The number of threads per block of the CUDA kernels depends on the maximum shared-memory size, number of registers per thread, and number of threads required to use full thread warps of the CUDA architecture. All these factors are directly connected to the complexity of the NLL evaluation, i.e. which PDFs are involved in the calculation and the dimension of the input data sample. Therefore this number depends on the user analyses. From our tests we have found a very small improvement (< 1%) of performance when we tune this number for each kernel. So, we have decided to simplify the procedure, using for all kernels a common value for the number of threads per block, independent of the tasks carried out by the kernels (the default value is 256). The number of blocks per kernel is then calculated from the number of events divided by the number of threads per block, rounded to the greatest integer number.
We would like to underline that data analysts can choose which implementation to use for the NLL evaluation among the original RooFit algorithm, the OpenMP or the CUDA implementations, by using a flag at runtime. They do not need to change their applications, i.e. there is a common interface for the user to the three implementations.

Tests
In the following tests we use a statistical model based on the BABAR analysis for the branching fraction measurement of the neutral B meson to η K 0 S decay [8]. The model has 3 variables and 5 species. In total there are 29 PDFs: 11 Gaussians, 5 polynomials, 3 Argus functions, which are combined in 5 PDFs for multiplication and for addition, respectively. In our fit we leave 16 parameters free to float. The events for each species are generated from the corresponding PDFs using Monte Carlo generation techniques, not considering the time for the generation in the tests (we only consider the time spent for the fitting procedure). We use RooFit v3.14, which is part of ROOT v5.28. As minimizer we use MINUIT2 of the same ROOT distribution [9].
The first test we perform is a comparison between the original RooFit algorithm and the OpenMP implementations, running with a single thread. We run the tests on a dual socket Intel Westmere-based system: CPU (L5640) @ 2.27GHz (12 physical cores, 24 hardware threads in total), 10x4096MB DDR3 memory @ 1333MHz. The system is running 64-bit Scientific Linux CERN 5.5 (SLC5), based on Red Hat Enterprise Linux 5 (server). The default SLC5 Linux kernel (2.6.18-194.8.1.el5) is used for all the measurements. We use the Intel C++ compiler version 11.1 (20100414). In this test we are also interested to see the effect of the vectorization in the new algorithm. To do that we compile the application deactivating the data auto-vectorization made by the compiler (flag -no-vec). The results are shown in table 1. Due to different optimizations applied by the compiler, we expect small differences on the results when running in the three different cases (well below the parameters error values), which are negligible from the perspective of the fit results (small fractions of the free parameter errors), but they can lead to a different number of NLL evaluations required for the minimization. For this reason we compare the results using the wall clock time divided by the number of NLL evaluations. We findin that the new algoritm implementation with a single thread is about 2.5x and 4.5x faster than the RooFit implementation when running without and with vectorization, respectively. The speed-up due to the data vectorization is 1.8x, close to the theoretical maximum of 2x, since the CPU can handle in parallel two floating point operations with double precision using SIMD instructions.
Then we run another test to see the scalability of the OpenMP implementation (with vectorization) when running in parallel with a fixed number of events (100,000). We use the same system described before. The threads are pinned to the cores running them. When we reach the maximum number of physical cores (12), then the threads are pinned to hardware threads (24), still maximizing the amount of physical cores used. Also in this case the number of NLL evaluations required for the minimization is different due to the parallel reduction approximation (ranging between 10336 and 12665). For this reason the scalability is calculated considering the time spent for each NLL evaluation. The results are shown in figure 1. We reach a 8.5x speed-up with 12 threads (maximum number of physical cores) and 9.9x with 24 threads (maximum load for the system). The main limitation in the scalability is given by the remaining sequential part of the program (mainly the calculation of the normalization integral of each PDF) and from handling in memory the multiple arrays of results.
The last test we present is the comparison between running the fit with OpenMP and CUDA implementations on the CPU and GPU, respectively. In this case we use the following system: single socket Intel Nehalem-based CPU (i7 965) @ 3.2GHz (4 physical cores), with 2048MB DDR3 memory @ 1333MHz, and Graphics Card ASUS ENGTX470 (based on the NVIDIA GF100 "Fermi" architecture). The operating system and the compiler versions are the same used in the previous tests, and we use CUDA v3.2. We do the test varying the number of events of the data sample. The comparison is done considering the time spent for each NLL evaluation, because of the different number of NLL evaluations required for minimization for the two different implementations. (the difference in number of calls is up to 12% of the total calls). We take as reference the optimized parallel OpenMP implementation, requiring four parallel threads, so that we fully load the available CPU. In case of the CUDA implementation, we include the time spent for the copy of the events from host memory to the device memory and for the copy of the array of final results back to the host memory. From the hardware point of view, we are comparing two systems which can be considered commodity systems: a single GPU, whose main target is for computer gaming, versus a standard single socket desktop system with 4 cores. The results are shown in figure 2. We can see how the CUDA implementation behaves better for high number of events, which is due to the specific ability of the GPU architectures to take advantage of multiple threads.

Conclusion
In this paper we have described a different algorithm for the NLL evaluation in maximum likelihood fits with respect to the algorithm used in the RooFit package. We implemented this algorithm to run in parallel on CPU, using OpenMP, and GPU, using CUDA. In our test the OpenMP implementation with a single thread is about 4.5x faster than the RooFit implementation (table 1). Furthermore the OpenMP algorithm was executed in parallel, giving a speed-up of about 10x with respect to a single thread execution in our test on 12 cores (24 hardware threads) system (figure 1). The comparison between the OpenMP and CUDA implementations are made using commodity systems, that can be considered, in terms of price and power consumption, easily accessible to general data analysts. In this case, running the OpenMP implementation in parallel (with 4 threads), we were able to reach a boost of 2.8x with the CUDA implementation (figure 2).