Parallel mining of time-faded heavy hitters

We present PFDCMSS, a novel message-passing based parallel algorithm for mining time-faded heavy hitters. The algorithm is a parallel version of the recently published FDCMSS sequential algorithm. We formally prove its correctness by showing that the underlying data structure, a sketch augmented with a Space Saving stream summary holding exactly two counters, is mergeable. Whilst mergeability of traditional sketches derives immediately from theory, we show that merging our augmented sketch is non trivial. Nonetheless, the resulting parallel algorithm is fast and simple to implement. To the best of our knowledge, PFDCMSS is the first parallel algorithm solving the problem of mining time-faded heavy hitters on message-passing parallel architectures. Extensive experimental results confirm that PFDCMSS retains the extreme accuracy and error bound provided by FDCMSS whilst providing excellent parallel scalability.

Introduction
In this paper we deal with the problem of mining in parallel time-faded heavy hitters (also called frequent items), and we present PFDCMSS, a novel message-passing based parallel algorithm which is a parallel version of the recently published FDCMSS sequential algorithm [4].
Mining of heavy hitters in a data stream has been thoroughly studied, and the problem is regarded as one of the most important in the streaming algorithms literature. Depending on the particular application, the problem is reported in the literature as hot list analysis [19], market basket analysis [2] and iceberg query [17], [1]. Even though there are many possible applications, we recall here some of the most important contexts to which the problem has been successfully applied: network traffic analysis [14], [16], [27], analysis of web logs [7], Computational and theoretical Linguistics [18].
All of the algorithms for detecting heavy hitters can be classified as being either counter or sketch based, the difference being that counter-based algorithms rely on a set of counters which are used to keep track of stream items, whilst sketch-based algorithms monitor the data stream by using a sketch data structure, often a bi-dimensional array data structure containing a counter in each cell. Stream items are mapped by hash functions to corresponding cells in the sketch. The former algorithms (counter-based) are deterministic, whilst the latter (sketch-based) are probabilistic.
Regarding counter-based algorithms, the first sequential algorithm has been designed by Misra and Gries [26]. Their algorithm was rediscovered, independently, about twenty years later by Demaine et al. [14] (this algorithm is known in the literature as the Frequent algorithm) and Karp et al. [22]. Among the developed counters-based algorithms we recall here Sticky Sampling and Lossy Counting [24], and Space Saving [25]. Sketch-based solutions include CountSketch [7], Group Test [12], Count-Min [11] and hCount [21].
Relevant parallel algorithms include [6], [3] and [5] which are message-passing based parallel versions of the Frequent and Space Saving algorithms. Shared-memory algorithms have been designed as well, including a parallel version of Frequent [31], a parallel version of Lossy Counting [30], and parallel versions of Space Saving [28] [13].
Recent shared-memory parallel algorithms for heavy hitters were recently proposed in [29]. Finally, accelerator based algorithms exploiting a GPU (Graphics Processing Unit) include [20] and [15]. Regarding related work, i.e., parallel algorithms specifically designed to solve the problem of mining time-faded heavy hitters, we are not aware of any other algorithm: to the best of our knowledge, ours is the first parallel algorithm solving the problem on messagepassing parallel architectures.
In this paper, we are concerned with the problem of detecting in parallel heavy hitters in a stream with the additional constraint that recent items must be weighted more than former items. The underlying assumption is that, in some applications, recent data is certainly more useful and valuable than older, stale data. Therefore, each item in the stream has an associated timestamp that will be used to determine its weight. In practice, instead of estimating items' frequencies, we are required to estimate items' decayed frequencies.
This paper is organized as follows. We recall in Section 2 preliminary definitions and concepts that will be used in the rest of the manuscript. We present in Section 3 our PFDCMSS algorithm and formally prove in Section 4 its correctness. Next, we provide extensive experimental results in Section 5, showing that PFDCMSS retains the extreme accuracy and error bound provided by the sequential FDCMSS whilst providing excellent parallel scalability. Finally, we draw our conclusions in Section 6.

Preliminary definitions
In this Section we introduce preliminary definitions and the notation used throughout the paper. We deal with an input data stream σ consisting of a sequence of n items drawn from a universe U; without loss of generality, let m be the number of distinct items in σ i.e., let U = {1, 2, . . . , m}, which we will also denote as [m]. Let f i be the frequency of the item i ∈ U (i.e., its number of occurrences in σ), and denote the frequency vector by f = ( f 1 , . . . , f m ). Moreover, let 0 < φ < 1 be a support threshold, 0 < < 1 a tolerance such that < φ and denote the 1-norm of f (which represents the total number of occurrences of all of the stream items) by ||f|| 1 .
In this paper, we are concerned with the problem of detecting in parallel frequent items in a stream with the additional constraint that recent items must be weighted more than former items. The underlying assumption is that, in some applications, recent data is certainly more useful and valuable than older, stale data. Therefore, each item in the stream has an associated timestamp that will be used to determine its weight. In practice, instead of estimating frequencies, we are required to estimate decayed frequencies. Two different models have been proposed in the literature: the sliding window and the time fading model. PFDCMSS works in the latter model. Furthermore, even though the basic ideas underlining the algorithm are also appropriate for an online distributed setting, here we are assuming that the entire dataset is available for offline processing.
The time fading model [23] [9] [8] does not use a window sliding over time; freshness of more recent items is instead emphasized by fading the frequency count of older items. This is achieved by computing the item's decayed frequency through the use of a decay function that assign greater weight to more recent occurrences of an item than to older ones: the older an occurrences is, the lower its decayed weight. Definition 1. Let w(t i , t) be a decayed function which computes the decayed weight at time t for the occurrence of item i arrived at time t i . A decayed function must satisfy the following properties: Related work has mostly exploited backward decay functions, in which the weight of an item is a function of its age, a, where the age at time t > t i is simply a = t − t i . In this case, w(t i , t) is given by , where h is a positive monotone non-increasing function. The term backward decay stems from the aim of measuring from the current time back to the item's timestamp. Prior algorithms and applications have been using backward exponential decay functions such as h(a) = e −λa , with λ > 0 as decaying factor.
In our algorithm, we use instead a forward decay function, defined as follows (see [10] for a detailed description of the forward decay approach). Under forward decay, the weight of an item is computed on the amount of time between the arrival of an item and a fixed point L, called the landmark time, which, by convention, is some time earlier than the timestamps of all of the items. The idea is to look forward in time from the landmark to see an item, instead of looking backward from the current time.
Definition 2. Given a positive monotone non-decreasing function g, and a landmark time L, the forward decayed weight of an item i with arrival time t i > L measured at time t ≥ t i is given by w(t i , t) = g(t i −L) g(t−L) . The denominator is used to normalize the decayed weight so that w(t i , t) is always less than or equal to 1 as requested by Definition 1.
Definition 3. The decayed frequency of an item v in the input stream σ, computed at time t, is given by the sum of the decayed weights of all the occurrences of v in σ: Definition 4. The decayed count at time t, C(t), of a stream σ of n items is the sum of the decayed weights of all the items occurring in the stream: The Approximate Time-Faded Heavy Hitters (ATFHH) problem is formally stated as follows.
Problem 1. Approximate Time-Faded Heavy Hitters. Given a stream σ of items with an associated timestamp, a threshold 0 < φ < 1 and a tolerance 0 < < 1 such that < φ, and letting g be a decaying function used to determine the decayed frequencies and t be the query time, return the set of items F, so that: • F contains all of the items v with decayed frequency at time t f v (t) > φC(t) (decayed frequent items); In the following, when clear from the context, the query time shall be considered an implicit parameter, so we write f v and C instead of f v (t) and C(t). The algorithm presented makes use of a Count-Min sketch data structure augmented by a Space Saving summary associated to each sketch cell. In the following, we recall the main properties of the Count-Min and the Space Saving algorithms in the case of non decaying frequencies, but the same properties also hold in a time-fading context.
Count-Min is based on a sketch whose dimensions are derived by the input parameters , the error, and δ, the probability of failure. In particular, for Count-Min d = ln 1/δ is the number of rows in the sketch and w = e/ is the number of columns. Every cell in the sketch is a counter, which is updated by hash functions. By using this data structure, the algorithm solves with probability greater than or equal to 1 -δ the frequency estimation problem for arbitrary items. The algorithm may also be extended to solve the approximate frequent items problem as well, by using an additional heap data structure which is updated each time a cell is updated. Since in Count-Min the frequencies stored in the cells overestimate the true frequencies, a point query for an arbitrary item simply inspects all of the d cells in which the item is mapped to by the corresponding hash functions and returns the minimum of those d counters.
Space Saving is a counter-based algorithms solving the heavy hitters problem. It makes use of a stream summary data structure composed by a given number of counters k n, n being the length of the stream. Each counter monitors an item in the stream and tracks its frequency. A substitution strategy is used when the algorithm processes an item not already monitored and all of the counters are occupied.
Let σ be the input stream and denote by S the summary data structure of k counters used by the Space Saving algorithm. Moreover, denote by |S| the sum of the counters in S, by f v the exact frequency of an item v and byf v its estimated frequency, letf min be the minimum frequency in S. If there exist at least one counter not monitoring any item,f min is zero.

The algorithm
In this section, we start by recalling our sequential algorithm FDCMSS [4]. The key data structure is an augmented Count-Min sketch D, whose dimensions d (rows) and w (columns) are derived by the input parameters , the error, and δ, the probability of failure. Whilst every cell in an ordinary CM sketch contains a counter used for frequency estimation, in our case a cell holds a Space Saving stream summary with exactly two counters. The idea behind the augmented sketch is to monitor the time-faded items that the sketch hash functions map to the corresponding cells by an instance of Space Saving with two counters, so that for a given cell we are able to determine a majority item candidate with regard to the sub-stream of items falling in that cell.
Indeed, by using a data structure S with two counters in each cell, and letting C i, j denote the total decayed count of the items falling in the cell D[i][ j], the majority item is, if it exists, the item whose decayed frequency is greater than C i, j 2 . The corresponding majority item candidate in the cell is the item monitored by the Space Saving counter whose estimated decayed frequency is maximum. We have proved that, with high probability, if a time-faded item is frequent, then, in at least one of the sketch cells where it is mapped, it is a majority item with regard to the sub-stream of items falling in the same cell. Therefore, our algorithm will detect it. Theorem 1. If an item i is frequent, then it appears as a majority item candidate in at least one of the d cells in which it falls, with probability greater than or equal to 1 − ( 1 2φw ) d . Regarding the error bound of our algorithm, let f i be the exact decayed frequency of item i in the stream σ and f i be the estimated decayed frequency of item i returned by FDCMSS. Let C be the total decayed count of all of the items in the stream. We have proved the following error bound.
Theorem 2. ∀u ∈ [m],f u estimates the exact decayed count f u of u at query time with error less than C and probability greater than 1 − δ.
The proofs of aforementioned theorems can be found in [4]. The algorithm's initialization requires as input parameters , the error; δ, the probability of failure; and φ, the support threshold. The initialization returns a sketch D. The procedure starts deriving d = ln 1/δ , the number of rows in the sketch and w = e 2 , the number of columns in the sketch. Then, for each of the d * w cells available in the sketch D we allocate a data structure S with two Space Saving counters c 1 and c 2 . Given a counter c j , j = 1, 2, we denote by c j .i and c j . f respectively the counter's item and its estimated decayed frequency. Finally, we set the support threshold to φ, select d pairwise independent hash functions h 1 , . . . , h d : [m] → [w], mapping m distinct items into w cells, and initialize the count variable, representing the total decayed count of all of the items in the stream, to zero.
Updating the sketch upon arrival of a stream item i with timestamp t i , shown in pseudo-code as Algorithm 1, requires computing x, which is the non normalized forward decayed weight of the item, and incrementing count by x. Then, we update the d cells in which the item is mapped to by the corresponding hash functions h j (x), j = 1, . . . , d by using the Space Saving item update procedure.
Let S denote the Space Saving stream summary data structure with two counters corresponding to the cell to be updated. Updating S upon arrival of an item works as follows. When processing an item which is already monitored by a counter, its estimated frequency is incremented by the non normalized weight x. When processing an item which is not already monitored by one of the available counters, there are two possibilities. If a counter is available, it will be in charge of monitoring the item, and its estimated frequency is set to the non normalized weight x. Otherwise, if all of the counters are already occupied (their frequencies are different from zero), the counter storing the item with minimum frequency is incremented by the non normalized weight x. Then, the monitored item is evicted from the Algorithm 1 Process Require: i, an item; t i , timestamp of item i; D, sketch data structure Ensure: update of sketch related to item i; update the local total decayed count. 1: procedure process(i, t i , D) 2: x ← g(t − t i ) compute the non normalized decayed weight of item i 3: lCount ← lCount + x update local total decayed count 4: for j = 1 to d do 5: SpaceSavingUpdate(S, i, x) update the sketch 7: end for 8: end procedure counter and replaced by the new item. This happens since an item which is not monitored can not have a frequency greater than the minimal frequency.
PFDCMSS, the parallel version of our sequential algorithm, works as follows. We assume the offline setting in which the stream items have been stored as a static dataset along with the corresponding timestamps. It is worth noting here immediately that our algorithm works in the streaming (online) setting as well. Indeed, in the former case (offline setting) we partition the input dataset and timestamps using a simple 1D block-based domain decomposition among the available p processes and then process in parallel the sub-streams assigned to the processes using Algorithm 1. In the latter case (online setting), we have instead p distributed sites, each handling a different stream σ i , i = 1, . . . , p processed again using Algorithm 1.
In the parallel version, once the sub-streams have been processed, one of the processes is in charge of determining the time-faded heavy hitters. In order to do so, all of the processes engage in a parallel reduction in which their sketches are merged into a global sketch which preserves all of the information stored in the local sketches. This sketch is then queried and the time-faded heavy hitters are returned.
In the distributed setting, one of the sites may act as a centralized coordinator or there can be another different site taking this responsibility. The coordinator broadcasts, when required, a "query" message to the p sites, which then temporarily stop processing their sub-streams, and engage in the sketch merge procedure. We can imagine the distributed sites as being multi-threaded processes, in which one thread executes Algorithm 1, temporarily stops when a query message is received from the coordinator, creates a copy of its local sketch and then resume stream processing whilst another thread engages in the distributed sketch merging procedure using the sketch copy.
In order to retrieve the time-faded heavy hitters, a query can be posed when needed. The query, shown in pseudocode as Algorithm 2, starts by determining the global decayed count for the whole stream σ. This requires a parallel reduction in which the local decayed counts are summed. It is worth noting here that the global decayed count is still non normalized; the normalization occurs dividing by g(t − L), where t is the query time and L denotes the landmark time. Then, we build, through a user's defined parallel reduction, a global sketch G which is obtained by merging the local sketches. To do so, each process invokes a parallel reduction by using the MergeSketch operator shown in pseudo-code as Algorithm 3.
The sketches are reduced as follows: for every corresponding cell in two sketches to be merged, the hosted Space Saving summaries are merged following the steps described in [5], i.e., building a temporary summary S C consisting of all of the items monitored by both S 1 and S 2 . To each item in S C is assigned a decayed frequency computed as follows: if an item is present in both S 1 and S 2 , its frequency is the sum of the its corresponding frequencies in each summary; if the item is present only in one of either S 1 or S 2 , its frequency is incremented by the minimum frequency of the other summary. At last, in order to derive the merged summary, we take only the 2 items in S C with the greatest frequencies and discard the others.
It is worth noting here that the sum of the counters in the stream summary data structure S related to a given cell D[i][ j] is equal to the value that the Count-Min sketch-based algorithm would store in the counter variable corresponding to that cell, i.e., the 1-norm of the frequency vector corresponding to the sub-stream falling in the cell through the pairwise independent hash functions. Thus, an augmented sketch is equivalent, from this perspective, to a Count-Min sketch and this property is preserved by the merge procedure. From now on we will call this property if c m . f g(t−L) > φ * gCount then 10: if p > φ * gCount then 12: return F 17: end procedure 1-norm equivalence.
However, merging Count-Min sketches simply requires adding the corresponding cells' counters. Indeed, via linearity, the sum of sketches is equal to the sketch of the sums. Instead, in our case, we need an ad hoc procedure in order to correctly merge the two Space Saving stream summaries hosted by the corresponding cells so that 1-norm equivalence property is preserved. Nonetheless, the augmented sketch which results from our parallel merge reduction is 1-norm equivalent to the Count-Min sketch obtained by summing the Count-Min sketches corresponding to our augmented sketches which are the input of the parallel merge reduction.
Once the global sketch is obtained, the query procedure initializes F, an empty set, and then it inspects each of the d * w cells in the sketch D. For a given cell, we determine c m , the counter in the data structure S with maximum decayed count. We normalize the decayed count stored in c m dividing by g(t − L), and then compare this quantity with the threshold given by φ * gCount (gCount being the normalized global decayed count). If the normalized decayed frequency is greater, we pose a point query for the item c m .i, shown in pseudo-code as Algorithm 4. If p, the returned value, is greater than the threshold φ * gCount, then we insert in F the pair (c m .i, p).
The point query for an item j returns its estimated decayed frequency. After initializing the answer variable to infinity, we inspect each of the d cells in which the item is mapped to by the corresponding hash functions, to determine the minimum decayed frequency of the item. In each cell, if the item is stored by one of the Space Saving counters, we set answer to the minimum between answer and the corresponding counter's decayed frequency. Otherwise (none of the two counters monitors the item j), we set answer to the minimum between answer and the minimum decayed frequency stored in the counters. Since the frequencies stored in all of the counters of the sketch are not normalized, we return the normalized frequency answer dividing by g(t − L).
At the end of the query procedure the set F is returned.

Correctness
Here, we prove that our algorithm correctly merges two FDCMSS sketches. The merge procedure preserves all of the properties of the sketch, including the fact that, considering the sum of the Space Saving counters in each sketch cell, an FDCMSS sketch is 1-norm equivalent to the classical Count-Min sketch.
It is worth noting here that we would obtain a correct result by using the merge procedure presented in [5] to combine the Space Saving summaries stored in the corresponding sketch cells, but we also want to impose 1-norm Algorithm 3 MergeSketch Require: D 1 , D 2 : sketches to be merged. Ensure: G, the merged sketch 1: procedure MergeSketch(D 1 , D 2 ) 2: for each S 1 i j ∈ D 1 , S 2 i j ∈ D 2 do 3: 15: end procedure equivalence, i.e., the additional condition that the sum of counters' values in each merged cell always reflects the total decayed count of the items which fell in the corresponding cells. Indeed, in [5] we showed how to merge Space Saving stream summaries in parallel. However, we have proved that our merge procedure satisfies the Space Saving properties described by eq. 2-5, and the following relaxed version of eq. 1: As shown in Theorem 3, which is the main result of this section, it turns out that k = 2 counters (i.e., majority item mining) is a special case: when the Space Saving summaries to be merged hold two counters, than the property in eq. 1 holds for the merged summary in its original form, that is |S| = ||f|| 1 , without modifying the merge procedure designed in [5]. Proof. The correctness of the parallel FDCMSS sketch merge algorithm derives from the correctness of the Space Saving merge procedure, already shown in [5]. It remains to show that, when looking to the sum of the Space Saving counters associated to each cell, the merged augmented sketch is still 1-norm equivalent to a Count-Min sketch, that is, the sum of the counters values is equal to the decayed count of all the items fallen in that cell.
Let us recall the merge algorithm for Space Saving summaries introduced in [5]. We will use the multiset notation, thus let us rewrite the properties of a Space Saving summary stated in equations 1-4, this time with reference to multisets. Indeed, we model the input stream as a multiset (also called a bag), which essentially is a set where the duplication of elements is allowed. We shall use a calligraphic capital letter to denote a multiset, and the corresponding capital Greek letter to denote its underlying set. In particular, we extend the traditional notion of multiset as follows. Instead of considering an indicator function which returns the multiplicity of an item, we use a function providing the decayed frequency of that item. Therefore, summing over all of the items we obtain the total decayed count in place of the cardinality of the multiset. N = (N, f N ) is a pair where N is some set, called the underlying set of elements, and f N : N → R is a function which provides the decayed frequency for each x ∈ N according to Definition 3.

Definition 5. A decayed multiset
The decayed count of N is expressed by whilst the cardinality of the underlying set N is From now on, when referring to either the exact or estimated frequency of an item, we shall mean the item's exact or estimated decayed frequency. Recall that our Space Saving stream summary data structure uses exactly k = 2 counters, and let N = (N, f N ) be the input decayed multiset, S = (Σ,f S ) the decayed multiset of all of the monitored items and their respective counters at the end of the sequential Space Saving algorithm's execution, i.e., the algorithm's summary data structure. Let |S| be the sum of the frequencies stored in the counters, f N (e) the exact frequency of an item e,f S (e) its estimated frequency andf S min the minimum frequency in S, wheref S min = 0 when |Σ| < 2. Indeed, even though a summary data structure has exactly 2 counters, it may monitor less than 2 items, since an item is actually monitored if and only if its counter's frequency is different from zero. The following relations hold, for each item e ∈ N: Now, let S 1 = (Σ 1 ,f S 1 ) and S 2 = (Σ 2 ,f S 2 ) be two summaries related respectively to the input sub-arrays N 1 = (N 1 , f N 1 ) and Theorem 3 in [5] states that if eqs. (10) -(12) hold for S 1 and S 2 and, if it is verified a relaxed version of eq. (9), i.e., it holds that then these properties continue to be true also for S M (it is worth noting here that eq. (13) also holds for summaries produced by the sequential Space Saving algorithm). The authors show that this is enough to guarantee the correctness of the merge operation, but, in general |S M | ≤ |N|.
In order to obtain S M , we start combining S 1 and S 2 to obtain S C , and then, if |Σ C | > 2, we take the two counters with the greatest frequency values in S C in order to build S M , otherwise we return S M = S C .
We can express the combine operation as shown by the following equation: In the special case of stream summaries holding exactly k = 2 counters, it holds that for i = 1, 2, |Σ i | ≤ 2, and |Σ C | ≤ 4. Now, suppose that |S i | = |N i |, (this is true when S i is produced by the sequential Space Saving) and let δ =f min S 1 +f min S 2 and x = |Σ C | − 2. Furthermore, suppose that the entries in S C are sorted in ascending order with regard to the counters frequencies.
As proved in [5], it holds that: where the sum is extended over the first x entries. We have to show that the difference xδ − x i=1f S C (e i ) is always equal to zero when k = 2, so that |S M | = |N|. When x ≤ 0, xδ = 0. In that case, S M = S C and |S C | = |N|. When x > 0, the first x counters of S C have values equal to δ. To see this, consider the two cases x = 1 and x = 2. When x = 2, that is, the two summaries to be merged contain different items and |Σ C | = 4, this is easily seen by simple computations: in fact, δ is the minimum value a counter in S C can assume, and there are at least two counters with this value in S C , obtained combining the two counters with minimum value in S 1 and S 2 . As a consequence these counters are the first two, and xδ − x i=1f S C (e i ) = 0. When x = 1, one of the following cases arises: 1. one of the summaries (without loss of generality, let us suppose it is S 1 ) contains two counters, the other summary (S 2 ) contains only one counter and no item is in common between the summaries. In this case, δ is equal to the minimum counter in S 1 sincef min S 2 = 0, but it is also the minimum counter in S C , hence it holds that δ −f S C (e 1 ) = 0 2. both summaries contain two counters and they have exactly an item in common. In this case we further have to distinguish three cases: (a) the item in common has the minimum frequency in both the summaries. The combined frequency of this item will be equal to δ which is the sum of the minimum frequencies of two summaries. Its combined frequency is also the minimum in S C , hence it holds that δ −f S C (e 1 ) = 0; (b) the item in common has the maximum frequency in both the summaries. Its combined frequency is also the maximum value in S C , and S C contains two distinct items with combined frequency equal to δ which is also the minimum in S C , hence δ −f S C (e 1 ) = 0; (c) the item in common appears with minimum frequency in one summary (without loss of generality, let us suppose in S 1 ) and with maximum frequency in the other summary (S 2 ). The combined frequency of the item which appears with minimum frequency in S 2 is equal to δ which again is the minimum frequency of the counters in S C , hence δ −f S C (e 1 ) = 0.
Taking into account that in all of the cases when x = 1 the S C contains at least one item whose combined frequency is equal to δ, it holds that xδ − x i=1f S C (e i ) = 0.
We have shown that all of the properties of a Space Saving summary of two counters are preserved by the merge procedure introduced in [5]. It suffices to guarantee that all of the properties stated for an FDCMSS sketch continue to hold after the parallel merge procedure depicted in the algorithm presented. In particular, it holds the property 1, which guarantees that a merged FDCMSS sketch continues to be 1-norm equivalent to a Count-Min sketch.

Experimental results
In this section, we report experimental results on synthetic datasets. Here, we thoroughly test our algorithm using an exponential decay function. All of the experiments have been carried out on the Galileo cluster machine kindly provided by CINECA in Italy. This machine is a linux CentOS 7.0 NeXtScale cluster with 516 compute nodes; each node is equipped with 2 2.40 GHz octa-core Intel Xeon CPUs E5-2630 v3, 128 GB RAM and 2 16 GB Intel Xeon Phi 7120P accelerators (available on 384 nodes only). High-Performance networking among the nodes is provided by Intel QDR (40Gb/s) Infiniband. All of the codes were compiled using the Intel C++ compiler v17.0.0.
Let f be the true frequency of an item andf the corresponding frequency reported by an algorithm, then the Relative Error is defined as ∆ f = | f −f | f , and the Average Relative Error is derived by averaging the Relative Errors over all of the measured frequencies.
Precision, a metric defined as the total number of true heavy hitters reported over the total number of candidate items, quantifies the number of false positives reported by an algorithm in the output stream summary. Recall is the total number of true heavy hitters reported over the number of true heavy hitters given by an exact algorithm. In all of the results we obtained 100% recall, even on a tiny sketch of size 4 x 800 (recall may be less than 100%, but this happens only when the sketch size is really minimal). For this reason, to avoid wasting space, we do not show here recall plots. Rather, we present Precision, Absolute Error, Average Relative Error (ARE), Updates/ms and runtime/performance plots since we are interested in understanding the error behavior and the algorithm's scalability when we use an increasing number of cores of execution. Table 1 reports the experiments carried out. For each different metric under examination, we varied n, the stream size in billions of items, ρ, the skew of the zipfian distribution, φ, the threshold and w, the number of sketch columns. All of the other parameters are fixed when varying one of the previous ones, and we show, on top of each plot, the fixed parameters' values.
Finally, we also present, for the metrics of interest, the results obtained by fixing the stream size and varying the number of cores utilized from 1 to 512. We conclude this section with a comparison between strong and weak scalability.
With the experiment 1 we aim at measuring the algorithm accuracy, the experiment 2 aims at measuring how the parallelization affects the algorithm's accuracy, finally experiment 3 is meant to measure the computational performance of the parallel algorithm measuring both strong and weak scalability. Table 1: Design of experiments. The input stream size n is expressed in billions, ρ denotes the skewness of the input data distribution, φ is the support threshold, w represents the number of columns in the sketch data structure, p is the number of cores and gs denotes the grain size, i.e. the number of elements (expressed in billions) in the input stream for each core. The nf value stands for not fixed

Algorithm accuracy
As shown by Figure 1, our parallel algorithm provides 100% Precision in all of the experiments carried out. Both the Absolute and the Average Relative Error, depicted respectively in Figure 2 and 3, have extremely low values, in particular with regard to their mean values. We observe that the Absolute Error is only slightly affected by the stream size n, and it's not affected at all by the threshold φ. The behaviour observed when varying ρ and w is expected. We observe a decrement of the Absolute Error in both cases, since, when the skew is higher, the number of frequent items in the corresponding zipfian distribution is lower and, when w is higher, increasing the sketch size provides better accuracy and, correspondingly, less error. Regarding the Average Relative Error, we observe the same qualitative behaviors in the experiments carried out.
Finally, the updates done per millisecond, shown in Figure 4, appear to be stable around 100,000 when varying the stream size n and the threshold φ. There is a visible increment (from 100,000 to 120,000) when varying ρ and a decrement (from 100,000 to 90,000) when varying w. These behaviors are expected for the same reasons we gave when analyzing the error. Indeed, processing a stream in which the number of frequent items is lower is usually faster. On the other hand, increasing the sketch size provides better accuracy, but more time is required to update the sketch.

Impact of the parallelization on the accuray
We now discuss the experimental scalability shown by our parallel algorithm. Figure 5 provides the results for the metrics under examination when testing strong scalability. That is, we fix the problem size (i.e., the stream size n) and increase the number of cores on which the algorithm is executed. As shown, Precision, Absolute and Average Relative Error are not affected at all by a strong scaling of the application, with Precision always equal to 100% and extremely low error values. Finally, the observed increment of the Updates/ms when varying the number of cores utilized is expected, due to the frequency updates made in parallel. Ideally, the throughput of the algorithm, measured as updates/ms, should increase with the same rate of increase of the number of MPI processes. Figure 6 is related to the comparison we did to test weak scalability, which refers to the scalability of a parallel application obtained when the problem size is increased along with the number of cores, so that we measure how the running time changes with regard to the number of cores for a fixed problem size per core (whilst, for strong scalability, we measure how the running time changes with regard to the number of cores for a fixed total problem size).

Computational perfomance
As shown, the plot for strong scaling is a log-log plot of the running time versus the number of cores. The dashed straight line with slope -1 indicates ideal scalability, whereas any upward curvature away from that line indicates limited scalability. The plot reports a good strong scalability even on 512 cores; this is due to the high number of items to be processed in the input stream which makes the computational time higher than the parallel overhead.
Regarding weak scalability, the corresponding plot provides an indication of loss of performance when scaling from 1 to 16 cores while we have a very good scalability from 16 cores up to 512. This can be explained considering the parallel architecture used for testing. One computing node is made of two octa-core Xeon processors, so that the cores share the main memory banks and the third level cache memory. In the weak scalability experiment the problem size increases linearly with the number of processes, hence also the total memory increases linearly; since from 1 to 16 cores we use only one computing node the memory contention between parallel processes increases. When varying from 16 to 512 cores we use several different computing nodes ranging from 1 to 32; each node runs 16 processes which compete for memory accesses as already discussed, hence the further slight loss of performance is due to the communication overhead.

Conclusions
We have presented PFDCMSS, a novel message-passing based parallel algorithm for mining time-faded heavy hitters, which, to the best of our knowledge, is the first parallel algorithm solving the problem on message-passing parallel architectures. We have formally proved its correctness by showing that the underlying data structure, is non trivially mergeable. However, the parallel algorithm is fast and simple to implement, and we have shown, through extensive experimental results, that PFDCMSS retains the extreme accuracy and error bound provided by FDCMSS whilst providing very good parallel scalability.


extfieldsofstudy:
Computer Science
Mathematics

text:
Subjective fairness: Fairness is in the eye of the beholder

We analyze different notions of fairness in decision making when the underlying model is not known with certainty. We argue that recent notions of fairness in machine learning need to be modified to incorporate uncertainties about model parameters. We introduce the notion of {\em subjective fairness} as a suitable candidate for fair Bayesian decision making rules, relate this definition with existing ones, and experimentally demonstrate the inherent accuracy-fairness tradeoff under this definition.

Introduction
Fairness is a desirable property of decision rules applied to a population of individuals. For example, college admissions should be decided on variables describing merit, but may also need to take into account the fact that certain communities are inherently disadvantaged. At the same time, individuals should not feel that another individual in a similar situation obtained an unfair advantage. All this must be taken into account while still caring about optimizing for a decision maker's utility function.
In particular, for a given distribution over a population, we wish to derive a decision rule that takes into account a merit variable, but also ensures fairness for members of disadvantaged groups. The problem becomes even more challenging when we take into account potential uncertainties in decision making models, which can even make strict notions of fairness impossible to satisfy.
As an example, consider the problem of fair prediction with disparate impact as defined by Chouldechova [2016]. Informally, their formulation defines a statistic a such that true category y (also called outcome or true label) is conditionally independent of a sensitive variable z given the statistic and the model parameters θ, i.e. y ⊥ ⊥ z | a, θ. When we face uncertainties in our modeling assumptions, the natural thing is to impose that the conditional independence holds if we marginalize the parameters out, i.e. y ⊥ ⊥ z | a. As we argue later in the paper, such a condition is impossible to satisfy, even if it holds for every possible parameter value, and we must incorporate subjectivity when model parameters are uncertain.
We instead develop a natural, and widely applicable framework for fairness that relies on the available information. We develop algorithms for achieving a few different notions of fairness within the subjective framework, and in particualr recently proposed concepts of fairness that are grounded in concepts of similarity and conditional independence. We argue that a suitable notion of similarity in the Bayesian setting is distributional similarity conditioned on the observations. For the latter, as independence is difficult to achieve uniformly in the Bayesian setting, we suggest a relaxation, for which we provide a small experimental demonstration. Chouldechova [2016] considers the problem of fair prediction with disparate impact. In this context, a statistic a : Ω → A is test-fair in some underlying probability space (Ω, Σ, P θ ), with respect to the outcome y : Ω → Y and sensitive variable z : Ω → Z, if y is independent of z under the statistic and parameter θ, i.e. if y ⊥ ⊥ z | a, θ. While the authors do not explicitly discuss the distribution P θ , it must be known in order for us to be able to find a statistic a satisfying this property. Kleinberg et al. [2016] consider feature vectors x : Ω → X and groups z : Ω → Z that correspond to the sensitive variables of Chouldechova [2016]. Their decision rule consists of a conditional distribution π(a | x) that assigns individuals to a finite A called "bins". There is also a score function U : A → [0, 1] corresponding to the utility of that bin. Finally, there is a categorical variable y : Ω → { 0, 1 }. They consider three fairness conditions, which we re-interpret below, subsuming two of them into a general "balance" condition:

Constraints on distributions and impossibility results
The authors show that these cannot can be simultaneously achieved under the distribution P π θ induced by the underlying latent parameter θ and the decision rule π. However, a sufficient condition for balance is the independence: U ⊥ ⊥ z | y, θ, π, i.e. when P π θ (U | y, z) = P π θ (U | y). In contrast to the condition in Chouldechova [2016], it is the decision variable U that is made independent given the ground-truth variable y, instead of the other way around. A stronger sufficient condition is the independence of the statistic a itself, rather than the utility variable.
Decision problems. While the aforementioned works focused on classification, Dwork et al. [2012] consider general decision rules π : X → A maximising expected utility, under a fairness constraint: where ∆ (P, Q) is a divergence between distributions P, Q, and ρ (x, y) is a metric in observation space. This formalizes a type of fairness where similar people are treated similarly. Even though there is no inherent notion of a sensitive variable, it's possible to capture it in the observation metric to a limited extent. In a similar decision-theoretic context, Corbett-Davies et al. [2017] consider tradeoffs between utility maximization and satisfaction of constraints on distributions, such as balance.
On sequential decision problems, such as multi-armed bandits and reinforcement learning,  define an algorithm as fair if it plays arms with highest means most of the time. Jabbari et al. [2016] study a similar notion for Markovian environments, whereby the algorithm's is fair the algorithm is more likely to play actions that have a higher utility under the optimal policy. However, this notion of fairness relies on oracle knowledge, while we focus on subjectivity. Here, y, z are clearly independent given a, θ. However, when we marginalize over θ they become dependent.

Fairness, subjectivity and learning.
Unfortunately it may may impossible to find a non-trivial decision rule under uncertainty, even for the case where only a single notion of fairness is appropriate. While such a notion may be satisfied for the correct underlying distribution P θ , it may not hold for an arbitrary estimator: In the Bayesian case, even if all P θ in a family satisfy a fairness condition, the marginal distribution may not.
Going back to our example from the introduction, we can easily see the impossibility of satisfying the disparate impact condition of Chouldechova [2016] from Figure 1. In this case, even if y ⊥ ⊥ z | a, θ, the same independence does not hold if θ is latent, no matter how a is distributed. This is clear from the graphical model, as y, z are connected through θ. When θ is observed, it is blocking and independence can be achieved.
In the Bayesian setting, even if all the members of our family satisfy a fairness condition, uncertainty about the correct distribution can create unfairness. In the following, we begin by a definition of the Bayes optimal rule, and then consider two achievable fairness constraints: smoothness and balance.

Fair Bayes decision rules
We begin with a simple statistical decision problem, where the decision maker observes x ∈ X , then takes a decision a ∈ A and obtains utility U (y, a) depending on an outcome y ∈ Y generated from some unknown distribution P θ (x | y).
We have some prior information β ∈ B in the form of a probability distribution on parameters θ ∈ Θ, an observation space X , an outcome space Y, and a family We consider a simple factorization, where β(θ | x, a) = β(θ), and a utility function U : Y × A → R that only depends on our action and the outcome. The decision diagram is given in Figure 2. Definition 1 (Bayes decision rule). The Bayes decision rule π : B × X → A is a deterministic policy that maximizes the utility in expectation, i.e. takes action where P β|x (y) = Θ P θ (y | x) dβ(θ) is the marginal distribution over outcomes conditional on the observations.
We can define analogues of the fairness concepts listed previously in terms of decision rules π, rather than any specific random variable. We will focus on Bayesian decision rules, i.e. rules whose decisions solely depend upon a posterior distribution.

Fairness as smoothness
The fairness notion in Dwork et al. [2012] requires action distributions to be similar for similar individuals, but the question of what "similar" is left open to the context. In the Bayesian setting, there is a natural candidate: the relevant information about an individual encoded in the distribution of outcomes P β|x (y):

Definition 2 (Smoothly fair Bayesian decision rule). A stochastic decision rule
where D (P Q) is some divergence between measures P and Q, and ∆ (P, Q) is a (possibly different) divergence. It is simply (k, β)-prior-fair with respect to the prior β, if The first definition ensures that, people whose merit distribution is similar under the same distribution β on parameters, will have a similar treatment. The second definition looks instead at how the treatment of one particular person who differ under a changing belief β on parameters.
We analyze Bayesian procedures that are commonly used in reinforcement learning problems. The first algorithm samples the underlying parameter θ from the posterior distribution and selects the decision that is optimal for the sample (i.e. marginalising over outcomes y). This is commonly known as Thompson samplingThompson [1933] in the reinforcement learning literatureWang et al. [2005], Agrawal and Goyal [2012].
Lemma 1. The Thompson sampling algorithm, i.e. making a decision a with probability equal to is prior-fair when ∆ is the KL divergence.
Proof. Thompson sampling can be seen as first drawing someθ ∼ β(θ | x) and then selecting the decision maximising Eθ(U | a) = Y U (y, a) dPθ(y | x), which depends only onθ. Through the post-processing inequality, we obtain the required result.
The second algorithm instead samples the latent outcome y, marginalising over parameters θ and corresponds to sampling actions according to their stochastic dominance.
Lemma 2. Stochastic dominance sampling, i.e. making a decision a with probability equal to is observation-fair when ∆ is the KL divergence.
Proof. This is in line with the proof of the previous lemma, but now we are drawing outcomes y instead of parameters θ.
For single-shot decisions, sampling has no informational advantage, as there is no exploration-exploitation trade-off. This makes the following result natural, but one should keep in mind that it also applies to fact that the (intractable) Bayes-optimal sequential decision rule is also deterministic.
Lemma 3. The Bayes-expected utility of stochastic dominance sampling is lower than that of Thompson sampling, which in turn is lower than that of the Bayes decision rule.
Proof. Let a * (x) max a E(U | a, x) for some arbitrary distribution. Then for any randomized policy π: . We first apply this to the distribution β(θ) to obtain the result for Thompson sampling. For stochastic dominance, note that P β (X > Y ) = Θ P θ (X > Y ) dβ(θ). As stochastic dominance can be implemented by first sampling a parameter and then sampling a dominant variable under this parameter, we can reapply this fact and obtain the final result.
We now consider a simple extension of Thompson sampling, which takes k parameter samples. As k → ∞, the algorithm becomes closer to the Bayesoptimal.
Proof. The first claim follows from the fact that we obtain k independent draws and from the additive property of the KL divergence. Optimality follows from standard concentration inequalities. For example, from Hoeffding's inequality, and bounded rewards, we obtain that with probability 1−δ our k-sample reward is ln 1/δ 2k -close to optimal. The result follows by choosing δ = k −1/3 , and using the boundedness of the utility.
Although these algorithms are defined for simple decision problems, they can be easily applied to sequential problems. This can be achieved if e.g. our decision space is the set of all behavioural policies. Now we look at another notion of fairness, which is instead places constraints on conditional independence of random variables.

Fairness as balance
The balance requirement in Kleinberg et al. [2016] involves equality in the expectation of a score function (depending on an observation x) under different values of a sensitive variable z, conditioned on the true (but latent) outcome y.
We can distinguish two cases for the dependency structure between x, y, z.
The first is a discriminative model satisfying y ⊥ ⊥ z | x, θ, our decision depends on x, but also on z, as x generates z. For a generative model satisfying x ⊥ ⊥ z | y, θ, we only have that x is not independent of z if θ is latent. However, it might still be possible to obtain a decision rule satisfying the following balance condition, analogously to Kleinberg et al. [2016]: Definition 3 (Balanced decision rule). A decision rule 1 π(a | x) is balanced with respect to P if a, z are independent under β for all y, i.e. if P π (a, z | y) = P π (a | y)P π (z | y), where P π is the distribution induced by P and the decision rule π.
Expanding the left hand side, we obtain P π (a, z | y) = P π (a | z, y)P (z | y) while for the right hand side, we have Equating the two terms, we obtain that This implies two sufficient conditions: The first is x ⊥ ⊥ z | y, but this is something that depends entirely on P , rather than the decision rule itself. Under the generative model of Figure 3b, this is achieved when θ is known, but not necessarily under β. The second is that the policy probabilities for every possible action a are orthogonal to the difference between the distribution of x for every possible choice of z, y.
Empirical formulation. For infinite X , it may be more efficient to rewrite (18) as This allows us to approximate the integral by sampling x, and can be useful for e.g. regression problems.

The optimally balanced decision rule.
Such a rule would maximize expected utility under the balance constraint. In our setting, however, we have no fixed P , but rather a belief β on a family of distributions P θ . This inherently makes our decision rule subjectively fair. We consider two variants of subjective fairness in this context: credible fairness and its natural relaxation, marginal fairness.
Credible fairness. In the Bayesian setting, we have a family of models { P θ } with a corresponding subjective distribution β(θ). As no single model is correct we need to consider fairness for every possible model. Rather than looking for a single useful decision rule π such that: ∀x, y, z, θ i.e. so that the constraint holds for all members { P θ | θ ∈ Θ } of the family, we instead choose to minimize some penalty function averaged over all models. For that reason, let us first define to be the deviation from balance for decision rule π under parameter θ. We now use this in the following definition: Definition 4. A decision rule is ǫ-credibly fair under β, and p, q ≥ 0, when ǫ ≥ Θ C(π, θ) dβ(θ).
In order to find a rule trading off utility for balance, we can maximize a convex combination of the expected utility and deviation. In particular, we can look for a rule π solving the following unconstrained maximization problem.
The rule maximising E π β f is a Bayes decision rule (Def. 1) for f under the distribution β. To find this, we can restrict ourselves to parametric decision rules, and use stochastic gradient descent. As in Section 3.1, if we only take k < ∞ samples to perform the gradient ascent over, the policy becomes an instance of k-Thompson sampling. When p = q = 2, the gradient procedure is quite simple, and is given in App. A, Marginal fairness. As an approximation, we can maximize with respect to the marginal model, i.e.
where E β , D β are defined in terms of the marginal distribution P β Θ P θ dβ(θ). This could be a useful formulation if calculating the marginal is easier than sampling.
A demonstration. Here we consider a discrete decision problem, for which we generate 100 observations. In particular |X | = |Y| = |Z| = |A| = 4, and U (y, a) = I {y = a}. The true model has a dependence of approximately 0.26. These observations are used to calculate a posterior distribution β(θ). The prior distribution is a simple Dirichlet-product, as the network is discrete, and we assume we know the structure, which is y → x → z. We then find two decision rules. The first uses the k = 10 samples from the full posterior distribution, and the second uses the marginal model. Figure 4 shows the estimated 2 and achieved 3 utility (SU, MU) and dependence (SD, MD) for the sampling-based and marginal decision rule respectively. Both rules are successful in trading off utility and dependency, with the rules based on the complete distribution being slightly better. The results match the estimated performance rather well for this toy problem.  Figure 4a shows the utility and dependence as measured during the optimization, while Figure 4b shows the actual utility and dependence for the marginal and credible decision rules for the true underlying distribution, as we vary λ.

Conclusion
Existing fairness criteria may be hard to satisfy or verify in a learning setting because they are defined for the true model. For that reason, we develop a subjective fairness framework, which can encompass existing criteria from a Bayesian viewpoint. This allows us to make some new connections, the stochastic gradient ascent algorithm, for example, can simultaneously satisfy similarity and independence definitions.
In addition, although the smoothness fairness criterion implies a stochastic decision rule, this rule can become deterministic when the information about individuals results in a singular distribution of the merit variable y. This implies a link between the independence and similarity conditions for certain information sets. In that case, it's also easy to see that the independence criteria will also be met. We believe that a further exploration of the informational aspects of fairness, and in particular in the Bayesian setting, will be extremely fruitful. Now ∇ w c w (y, z) ⊤ c w (y, z) = ∇ w a c w (y, z) 2 a (26) = a 2c w (y, z) a ∇ w c w (y, z) a ∇ w c w (y, z) a = x ∇ w π w (a | x)D(x, y, z), while Combining the two terms, we have − 2λ y,z c w (y, z) a D(x, y, z) dΛ(x), .
where Λ is the Lebesgue measure. We now derive the gradient for the ∇ w π w term. We consider two parameterizations.
Independent policy parameters. When π(a | x) = w ax , we obtain ∂π(a ′ | x ′ )/∂ax = I {ax = a ′ x ′ }. This unfortunately requires projecting the policy parameters back to the simplex. For this reason, it might be better to use a parameterization that allows unconstrained optimization.
Softmax policy parameters. When π(a | x) = e wax / a ′ e w a ′ x .


extfieldsofstudy:
Mathematics
Physics
Computer Science

text:
Minimization of Antenna Quality Factor

The optimal currents on arbitrarily shaped radiators with respect to the minimum quality factor Q are found using a simple and efficient procedure. The solution starts with a reformulation of the problem of minimizing quality factor Q as an alternative, so-called dual, problem. Taking advantage of modal decomposition and group theory, it is shown that the dual problem can easily be solved and always results in minimal quality factor Q. Moreover, the optimization procedure is generalized to minimize quality factor Q for embedded antennas, with respect to the arbitrarily weighted radiation patterns, or with prescribed magnitude of the electric and magnetic near-fields. The obtained numerical results are compatible with previous results based on composition of modal currents, convex optimization, and quasi-static approximations; however, using the methodology in this paper, the class of solvable problems is significantly extended.

I. INTRODUCTION
T HE quality factor Q (or Q-factor) is of primary importance in the design and analysis of electrically small antennas because of its approximate proportionality to fractional bandwidth [1]- [3]. Classical approaches to estimate the value of an antenna's quality factor Q center around evaluating the frequency sensitivity of the antenna in different ways: from the antenna's electric field integral equation [4], equivalently from its impedance (reactance) matrix [5], or directly from its input impedance [6], [7]. In many cases, these techniques can reasonably predict the bandwidth potential for high-Q antennas [1], [7], [8].
The determination of the fundamental bounds on antenna performance has evolved from methods studying canonical bodies, such as spherical shells [9], to the more general analyses of arbitrarily shaped radiators [3]. Much of this work is motivated by seeking an optimal radiator's performance quantified in terms of quality factor Q, antenna gain G, and radiation efficiency η rad as functions of electrical size [1], Manuscript  M. Gustafsson is with the Department of Electrical and Information Technology, Lund University, 221 00 Lund, Sweden (e-mail: mats.gustafsson@eit.lth.se).
K. Schab is with the Department of Electrical and Computer Engineering, Antennas and Electromagnetics Laboratory, North Carolina State University, Raleigh, NC, USA (e-mail: krschab@ncsu.edu). [2]. Certain antenna parameters (or their combinations) can be analyzed in the static limit [10] or be reformulated as convex optimization problems for larger structures [11], [12]. Effects of finite ground planes, lossy materials, or magnetic currents can also be included [5], [12], [13]. However, the unconstrained minimization of quality factor Q has not been formulated as a convex problem, making its solution much more difficult.
The method employed in this paper proposes a dual problem formulation or, in other words, a relaxation of the unconstrained optimization problem [12]. Determining the optimal perturbation parameter while taking the smallest eigenvalue of the related generalized eigenvalue problem forms a problem which is directly solvable via the bisection algorithm [14]. Because the original problem of minimizing quality factor Q is generally non-convex, the difference between the dual solution and the solution to the original problem (i.e., the duality gap) can be non-zero [11], [12]. However, we show in this paper that through proper treatment of degenerate eigenmodes within the dual problem the duality gap can always be made zero. Taking advantage of the von Neumann-Wigner theorem [15] and work on crossing avoidance [16] and group theory [17], [18] we also determine when such degeneracies occur based on symmetry.
The presented approach is generalized to several related but previously unsolved antenna problems. The minimum quality factor Q for antennas near finite ground planes is found by dividing the antenna body into controllable and uncontrollable sub-regions [5]. It is also shown that minimization of the stored energy with arbitrary weighted constraints on the radiation intensity and the power in the electric and magnetic near-fields can be solved in similar fashion. This gives physical bounds for several antenna problems. Particularly, usage of the radiation intensity in a single direction generalizes the bounds on the partial gain in [3], [5] to the gain.
The proposed minimization procedure is computationally inexpensive and easy to implement. The results are also comparable with recently published methods using modal decomposition [19]- [21] to determine the optimal currents with respect to the minimum quality factor Q. For simplicity, the presented analysis is restricted to compactly supported electric surface current densities in free space. The results can directly be generalized to volumetric and magnetic current densities [22], [23].
The paper is organized as follows. The optimization procedure used to determine the minimum quality factor Q is introduced in Section II and its important properties are discussed in Section III, including a proof of method's convergence to the global minimum. The method is generalized in Section IV towards the utilization of advanced constraints. Validity of the method is verified by examples of various complexity in arXiv:1612.07676v1 [physics.class-ph] 22 Dec 2016 Section V. The consequences of the method are discussed in Section VI and the paper is concluded in Section VII.

II. MINIMIZATION OF QUALITY FACTOR Q
In this section, the problem of minimizing quality factor Q is formulated, along with the associated dual problem and its solution.

A. Problem Statement
For a given surface or volume of compact support Ω, we would like to determine the minimal quality factor Q lb , defined as Q lb = min in which I is yet unknown vector containing expansion coefficients of the current density J (r), r ∈ Ω, represented in a suitable basis {ψ n } [24], [25], i.e., The quality factor Q in (1) is defined as [26] Q For a given current distribution represented by I, the stored magnetic and electric energies are approximated as [4], [5] W m ≈ 1 8 respectively, and the radiated power P r is given by Here, Z = R + j (X m − X e ) is the impedance matrix based on the Galerkin method [24], X = X m − X e , ω is angular frequency, and superscript H denotes the complex conjugate transpose. The quadratic forms (4a) and (4b) based on the stored energy matrices X m and X e accurately predict the stored energies W m and W e for currents with support up to the electrical size of at least ka = 1, in which k is the wavenumber and a is the radius of smallest sphere circumscribing Ω [12]. These matrices can easily be obtained from method-of-moments codes [27].

B. Formulation and Solution of Dual Problem
Following [12], the original problem in (1) can be relaxed by convex combination 1 of Q m and Q e from (3), i.e., for any current I, for 0 ≤ ν ≤ 1. Minimizing both sides of (6) over I gives the lower bound Q lb in (1) on the LHS, and the solution to the dual problem for a given value of the perturbation parameter, ν, on the RHS of (6), i.e., We denote the solution to the RHS of (7) as Q ν and which can be easily obtained by finding the smallest eigenvalue solution to Inserting the solution Q ν into the RHS of (7) and recognizing that Q (I ν ) from (3) must be greater than or equal to the lower bound Q lb , we arrive at the three-part inequality where the following notation is adopted: Using these quantities, Q ν is also written as a convex combination If equality holds between the far-left and far-right sides of (10), then the duality gap is zero and Q lb is effectively found by solution of the dual problem. In the following section, we show that this is always the case. Hence, the concave dual problem of maximizing Q ν over ν can always be used to find the solution to the unconstrained minimization problem, i.e., A straightforward approach to solving the above problem is to use the bisection algorithm [14], where for each evaluated ν the value of Q ν is obtained by finding the smallest eigenvalue solution to (9).

III. PROOF OF A TIGHT LOWER BOUND
The following discussion examines how the properties of the eigenspace associated with I ν and Q ν always lead to equality in (10) for some value of ν.

A. One-dimensional Eigenspace
If the quality factor Q ν has one continuous associated eigenvector for all ν and at some ν 0 this eigenvector is self resonant, i.e., Q mν0 = Q eν0 , then equality holds in (10) at ν = ν 0 . Figure 1a depicts this self-resonant point for an L-shaped plate with largest dimension L = 0.1λ, where λ denotes the wavelength. By (10), Q lb is given by the value of Q eν and Q mν at their intersection, shown as a red dot in Fig. 1a.
If the quality factor Q ν has one continuous associated eigenvector for all ν but self-resonance is never achieved, (implying Q mν > Q eν for all ν or Q eν > Q mν for all ν) then equality in (10) must hold at either ν = 0 or ν = 1. Hence, when Q ν is always associated with a one-dimensional eigenspace, there exists a ν which leads to equality in (10).
The electric and magnetic modal quality factors Qmν and Qeν from (11) are depicted for two different shapes. The L-shape without any degeneracy is depicted on the left part. As all the curves are continuous, the minimal quality factor Q lb can be found for such ν where Qmν = Qeν . The rectangular plate is depicted on the right part of the figure. The minimal quality factor Q lb can be found in case the degenerate eigenspace is taken into account.

B. Instantaneous Degeneracies (Crossings)
In cases where Q ν has two or more degenerate eigenvectors at a discrete value ν 0 , the right hand side of (10) remains single-valued while the left hand side of (10) is not uniquely defined. This is because at the degeneracy I ν can be any arbitrary vector within the multidimensional eigenspace associated with Q ν . We can always choose a vector from this eigenspace which satisfies the self-resonant condition, provided the eigenvectors at ν 0 − δν and ν 0 + δν have opposite ordering of Q e and Q m . Satisfying the self-resonance condition in this way implies equality in (10) at ν 0 . In the specific case of a twodimensional eigenspace with basis vectors I 1 and I 2 , the selfresonant eigenvector I ν,sr is given by where the real-valued coefficient χ is the solution to and with the superscript T denoting the matrix transposition. Note that discrete degeneracies such as these occur only when a surface's geometric symmetries are precisely represented [15], [17], [18], in which case ∆ 12 , by way of symmetry, is zero if I 1 and I 2 are chosen to be continuous with the eigenvectors at ν 0 − δν and ν 0 + δν. From this we see that the phase associated with the combination in (14) is arbitrary.
If the eigenvectors at ν 0 − δν and ν 0 + δν do not have opposite ordering of Q e and Q m , then (15) has no real valued solution and the degeneracy can be otherwise ignored. In this case, the rules for Q ν with one-dimensional eigenspaces can be used to determine that there exists some ν where equality holds in (10). (17) Notation Bilinear form Ref.

Mode expansions
Ohmic losses Pσ Tab. I. Summary of additional quadratic forms which can be used to further constrain the radiated field, lost power or the near-field. Figure 1b shows Q mν and Q eν for a 0.1λ × 0.05λ rectangular plate. Due to the geometric symmetry of the plate, Q ν has a multidimensional eigenspace at ν ≈ 0.82, where Q mν and Q eν are discontinuous. At this value of ν, we can apply the preceding method of finding a self-resonant I ν within the degenerate eigenspace. Using this self-resonant vector, the value of Q lb is given by (10), shown as a red dot in Fig. 1b.

C. Distributed Degeneracies
If Q ν has a multi-dimensional eigenspace for all ν, then I ν should be selected as a vector which minimizes max{Q mν , Q eν } within the eigenspace. With the eigenvector defined in this way, the rules for Q ν with one-dimensional eigenspaces can again be used to determine that equality must hold in (10) for some ν. This case is very rare in practice and will only occur when certain geometric symmetries are precisely represented, e.g., a square with rectangular mesh or a sphere with global basis functions. The latter of these is discussed further in Section V-C.

IV. GENERALIZATION OF THE MINIMIZATION PROBLEM
The results of the previous sections are valid for arbitrary positive semi-definite matrices X m , X e , R and are hence directly generalized to several important antenna problems. For notational convenience, we rewrite the generalized eigenvalue problem (9) as In addition to the radiated power and stored energy, we consider quadratic forms enumerated in Table IV. These quantities are combined with the stored energy (4a) and (4b) and radiated power (5) to form eigenvalue problems for antennas with prescribed radiation patterns, maximization of G/Q, near-field constraints, and tradeoffs between stored energy and ohmic losses.
The classical bounds on the quality factor Q are formulated for the electric dipole (TM), magnetic dipole (TE), and combined cases [1]- [3], [9]. The TM and TE cases are determined from an expansion of the radiated field in spherical modes and division of the radiated power (5) in its TM and TE parts. Expansion of the radiated fields in spherical modes defines row matrices F α with elements [5] where Z 0 is the free space impedance, k = 2π/λ is the wavenumber, v α denotes the real valued regular spherical vector waves [28], α is the multi-index of the mode and ψ n are the basis functions used in (2). The decomposition of TE and TM modes corresponds to a decomposition of the resistance matrix in its TE and TM parts R = R TE + R TM . The eigenvalue problem for minimizing quality factor Q for antennas radiating only TM fields is (17) with The radiation pattern and the radiated power can also be expressed in the partial radiation intensity U (r,ê) for the directionr and polarizationê. The partial radiation intensity U (r,ê) is further rewritten using the far-field component e * · F (r) = F(r,ê)I, where the complex conjugate is used for elliptic and circular polarized cases. The row matrix F has the elements and The minimum stored energy for an antenna having the radiation intensity weighted with w 1 and w 2 is given by (17) with (22) whereê 2 =r ×ê * 1 . An alternative interpretation for the case with w 1 (r) and w 2 (r) as the characteristic functions for an angular region is as the minimum quality factor Q for antennas with negligible radiation outside the angular region.
The particular choice of w 1 (r ) = w 2 (r ) = δ(r −r) in (22) gives an optimization problem for the quotient between quality factor Q and antenna gain G The optimization problem for lower bounds on Q/G or equivalently upper bounds on G/Q is hence reformulated as minimization of where the spherical unit vectorsθ andφ are used to express the polarization. The generalized eigenvalue problem (17) contains the rank two matrix It is also possible to use the dissipated power in the near field. Near fields are of interest for exposure, wireless power transfer, mutual coupling, and isolation [29]. Here, the matrix R is expressed as quadratic forms constructed by the nearfield matrices N e,n and N m,n that are used to expressed the electric and magnetic near fields as E(r n ) = N e,n I and H(r n ) = N m,n I, respectively.
Ohmic losses can be included to investigate bounds on the radiation efficiency η rad of antennas. For lossy antennas it desired to have both low losses and low stored energy for a fixed radiated power. This leads to a multi-objective optimization problem which can be analyzed using Pareto fronts [30] by minimization of for fixed value of the parameter 0 ≤ ς ≤ 1, where I H R σ I is constructed to represent relevant dissipative losses in the system [21], [25]. This problem has the same form as (17) with the substitutions R = R and and is hence easily solved. The parameter ς determines the weight between the stored energy and losses. Small and large values of ς emphasizes low losses, expressed as η −1 rad − 1, and low Q/η rad , respectively. There are several other interesting multi-objective formulations, e.g., combining the Q/G in (24) with 1/G in (23) as can be used to investigate the trade-off between quality factor Q and G. This is similar to the convex optimization formulation for minimum quality factor Q for superdirective antennas [5].

V. RESULTS
In this section, the procedure to calculate the minimum quality factor Q lb is demonstrated on a series of numerical examples of various complexity using bodies made of perfect electric conductor (PEC). Important consequences are discussed in Section VI.
First, a structure without geometric symmetry is treated in Section V-A. It is shown that these cases are easily solvable, cf. Fig. 1a. Second, a structure possessing geometric symmetries is analyzed in Section V-B and the solution is found in the degenerate eigenspace, cf. Fig. 1b. The spherical shell is considered in Section V-C as the third example, being a representative of canonical bodies which are solvable analytically [20]. The spherical shell is of great interest here as it possesses the highest degree of symmetry in R 3 . The last example in Section V-D deals with a simple structure subject to sub-domain optimization.
All dimensions are normalized to kL or, whenever proper, to ka, in which L is a length of the structure (to be specified for each example individually) and a is the radius of the smallest sphere completely surrounding the sources.

A. L-shaped Plate
A planar, L-shaped PEC structure ( ) of length L, width L/2, cut of L/2 × L/4 and electrical size kL ≈ 0.628 is discretized into 792 triangular elements and considered here as the first example, see the inset in Fig. 2.
Because the structure possesses no geometrical symmetries, the procedure described in Section II can be followed without special treatment for degenerate modes. In order to find the parameter ν at which Q lb = Q lb , the bisection algorithm [14] can be employed to solve 2 the maximization problem in (13). The minimum quality factor Q lb occurs at ν ≈ 0.865, see Fig. 2. It reaches Q lb = 128 and yields Q mν = Q eν , see Fig. 3 as well. The radiation patterns of the current I ν are shown in Fig. 3 for ν = 0.75 and ν = 0.95, i.e., for values where Q mν < Q eν and Q mν > Q eν , respectively.
The behavior of the separate Q mν and Q eν terms is depicted in Fig. 3. It can be seen that in the vicinity of the optimal solution the quality factors Q mν and Q eν change abruptly but still smoothly. The sharp minimum of quality factor Q lb in Fig. 2 is direct consequence of this rapid change. Notice that because no symmetries are present on this structure, the same results can be obtained by, e.g., employing rectangular mesh elements.

B. Rectangular Plate
A structure with geometrical symmetries, a rectangular plate, is treated as the second example. This object remains unchanged under two axes of reflection and 180 o rotation. Hence, it has a non-trivial symmetry group [18] and, as such, a degenerate eigen-space may occur during the solution of (13).
The PEC plate ( ) has dimensions L × L/2 and electrical size kL ≈ 0.628. It is divided into 1540 rectangular elements and calculated in a method-of-moments code tailored to preserve the geometric symmetry of the plate.
When solving (13), a degenerate eigenspace is found at the maximum of Q ν , see Fig. 4. As a consequence of the degeneracy, it is not immediately obvious that Q lb = Q lb at this point. However, treating the degenerate eigenspace according to (14)- (16) one finds that this equality holds when the self-resonant eigenvector within the eigenspace is considered; leading to Q lb = Q lb = 103. The existence of the degenerate eigenspace introduces the additional degree of freedom for determination of the optimal current. This additional parameter, denoted in (14) as φ, is the phase between degenerate eigenmodes and its importance is studied analytically in the next example. The appearance of the degeneracy in this example results from the strict enforcement of the plate's symmetry. By perturbing the simulation through rounding errors or by using triangular basis functions [25], [31] which may not maintain symmetry, the degeneracy can be effectively eliminated.

C. Spherical Shell
The PEC spherical shell of radius a is examined in this section. Spherical harmonics [32] are used as basis functions for which the impedance matrix and the X m , X e , and R matrices are diagonal and can be evaluated analytically [33].
The maximization procedure (13) is recapitulated in Fig. 5 for the particular choice of ka = 0.5. Since the eigenspace (9) of the spherical shell is degenerated, a gap of width Q lb − Q lb occurs. It can be, however, effectively removed by using procedure from Section III. Interestingly, Q mν and Q eν are constant with respect to ν, which means that Q ν is a linear function of ν. This is verified in Fig. 5 and simplifies locating the minimum quality factor Q lb ≈ 9.73, cf. [20], [34]. The second case depicted in Fig. 5, deals with maximization of G/Q, here reformulated as minimization of Q/G. There is no duality gap for this case, in agreement with the partial G/Q case [12] and the procedure gives Q/G ≈ 3.31 (G/Q ≈ 0.302).
The minimal quality factors Q lb and Q G/Q are depicted in Fig. 6 as functions of the sphere's electrical size, ka. The quality factor Q G/Q is calculated from (3) with the current minimizing the ratio G/Q. The quality factors are depicted for the cases of arbitrary radiated fields (TM+TE), radiated electric modes (TM), and radiated magnetic modes (TE) calculated using (19). The computed quality factors agree with the results in [20]. We observe that Q G/Q ≈ Q lb for ka 1 and the slight increase of Q G/Q over Q lb for larger ka is due to excitation of higher order spherical modes, i.e., Q G/Q = Q lb if the radiated fields are restricted to the dipole  Fig. 6. Minimal quality factors Q for the spherical shell of radius a studied in dependence on electrical size ka. The blue curve depict minimal quality factor Q lb , the orange curve depict quality factor Q Q/G , corresponding to the minimal value of Q/G. All curves are also categorized with respect to different radiation constrains: TM+TE stands for combination of electric and magnetic radiation modes, TM stand for electric (dipole-like) and TE for magnetic (loop-like) modes. All quality factors are normalized to (ka) 3 in order to flatten the curves.
modes. Moreover, as the phase between modes can be chosen arbitrarily for the degenerate eigenspace (14) the G/Q ratio can be increased by proper selection of TM-TE mode doublet with optimal phase shift. The spherical shell is classical example of a symmetrical body possessing distributed degeneracies. As explained in Section III-C, this gives us additional freedom selecting the preferable optimal currents. For non-symmetrical body there is always unique current yielding minimum quality factor Q and having typically elliptical polarization. However, considering body with distributed degeneracies, we can freely select preferable polarization.

D. Current Optimization in Sub-domain Region
The last example considers sub-domain antenna optimization [5], where the region Ω is divided into a controllable antenna region Ω A and an uncontrollable ground region Ω G . Currents in Ω A (denoted I A ) may be controlled via excitation voltages, whereas currents Ω G (denoted I G ) are restricted to those induced by currents in Ω A . We have I = I A ∪ I G and Ω = Ω A ∪ Ω G , and consequently the discretized EFIE can be written in block matrix form as leading to Leveraging the division of controllable and uncontrollable currents and the relation in (30), any quadratic form of the total current can be compressed into a quadratic form of the controllable current, i.e., where for any matrix A operating on the current vector I.
To investigate sub-domain optimization thoroughly, several different controllable regions are considered, see Fig. 7 and Table II. In all cases, the mesh refinement of the controllable region Ω A is improved to provide enough degrees of freedom for optimal current. The RWG basis functions coinciding with the boundary between Ω A and Ω G were considered as controllable.
The rectangular plate of the same dimensions as in Section V-B is treated in Fig. 7 and the dependence on size of controllable region and electrical size is studied. It can be seen that the dependence on the electrical size embodies the same course for different areas of region Ω A . Practically the same quality factors Q lb for curves c/L = 0.1 and d/L = 0.2 in Fig. 7 suggest that there is no significant difference between various sub-domain divisions while the constant ratio between areas Ω A and Ω G is preserved. It is, however, demonstrated in Table II that this observation is not completely trueinterestingly, if the controllable region Ω A is well-chosen, the loss of control over 50% of the total area led to only 30% increase of minimal quality factor Q lb . Of particular interest is the third candidate in Table II, for which the controllable region covers the outer loop of the rectangle. This region seems to be the optimal choice as we know that it resembles the loop-like optimal current, cf. Fig. 4. However, the value of the quality factor Q lb is quite poor due to capacitive coupling between controllable and induced current. This is important observation, especially in the context of small antenna design, since often the parasitic ground plane is prescribed and the controllable region for antenna design is limited. Results presented in this work imply that there exists both an optimal shape and optimal placement of the controllable area.

VI. DISCUSSION
It has been proved in Section III and verified on practical examples in Section V that solving the maximization problem in (13) always minimizes quality factor Q. In Section IV, other problems such as minimization of the Q/G ratio, were discussed. Several points not yet discussed regarding the implementation and properties of this method are briefly treated in the remaining part of the paper.

A. Requirements
The impedance matrix Z must be symmetric, i.e., Z = Z + Z T /2. This is always true if Galerkin testing is adopted and properly implemented [25]. Fortunately, this is a common choice when formulating the method of moments. In addition, it is required that the real part of the impedance matrix is positive semi-definite, i.e., R 0. This last requirement is grounded in the physical argument that any current distribution must radiate non-negative power. Due to numerical errors, many method of moments codes produces a matrix R with small negative eigenvalues. These negative eigenvectors should be removed to enforce R 0, [12].
From all presented results, we observe that the optimal current which minimizes quality factor Q tends to equalize magnetic and electric energies, i.e., Q eν = Q mν . Whether explicitly mixed at a discrete degeneracy or implicitly occurring at a smooth maximum when solving (13), we can interpret the balance of magnetic and electric energies as blending of predominately electric (dipole-like) and magnetic (looplike) type currents. It is important to note that the choice of discretization can directly impact the ability to represent looptype currents, e.g., a straight wire dipole represented with onedimensional basis functions (no circumferential flow) lacks the ability to support divergence-free loop currents.  Table II for the same electrical size (kL ≈ 0.628). The optimal currents are real-valued. The arrows in the quiver plots have normalized length and they depict the direction of the current. The color maps depict the normalized absolute value of the current density at each triangle. The density of mesh grid was scaled down in order to produce lucid illustrations.

B. Properties
As compared to other methods, it can be somewhat surprising that a purely real valued current is sufficient to deliver the minimum quality factor Q, cf. [21]. This is, however, well-aligned with recent attempts utilizing, e.g., characteristic modes, [20]. With respect to those modal approaches, this method can be understood in the same way -as a method combining modes. However, because typically only one or two modes are needed, we can conclude that the solution here is represented in its natural basis.
The minimization procedure can be applied for any problem with the convex combination of stored energy operators on the LHS and a positive semi-definite matrix on the RHS of [12]. Furthermore, and contrary to the previous attempts, the weighted matrix on the RHS is not restricted to being a rankone matrix, therefore the class of antenna problems solvable by convex optimization is greatly enlarged. Notice, however, that the stored energy operators X m and X e in (4a) and (4b) fail to deliver the true stored energy for higher frequencies [35], i.e., ka > 1. In the context of this work, this fact is only a minor technicality and does not restrict the generality or applicability of the proposed method.
Finally, notice direct portability of the method to magnetic current or combination of electric and magnetic currents.

VII. CONCLUSION
The long-persistent problem of quality factor Q minimization for arbitrarily shaped body has been reformulated in convex form. Consequently, it has been solved using straightforward and computationally inexpensive procedure. As long as the required algebraic properties are guaranteed, the same procedure is directly applicable to other problems as well. Namely, the Q/G ratio can easily be maximized (rewritten as minimization of Q/G), problems involving various near-or far-field restrictions can be solved, ohmic losses can be taken into account and the presence of regions with uncontrollable currents are now solvable. The enumerated problems are solved by the same formalism, which makes the potential implementation appealing for all scientists and other workers in the field.
One of the byproducts of the method is the link established between several diverse techniques which have been merged and further generalized. An important outcome of this work is the understanding that the current optimal with respect to quality factor Q and G/Q have typically a uniquely defined inter-modal phase. For the special case of symmetric objects, this phase may be arbitrary for some quantities being minimized.
The presented approach also opens new questions to be addressed. For example, it seems clear that, when having control of currents over only part of a structure, there is an optimal shape of the controllable area to have the greatest impact on the antenna's performance. Another example is the open problem of controlling interactions within multiantenna systems and the existence of optimal currents for these scenarios.