Bit-Vector Model Counting using Statistical Estimation

Approximate model counting for bit-vector SMT formulas (generalizing \#SAT) has many applications such as probabilistic inference and quantitative information-flow security, but it is computationally difficult. Adding random parity constraints (XOR streamlining) and then checking satisfiability is an effective approximation technique, but it requires a prior hypothesis about the model count to produce useful results. We propose an approach inspired by statistical estimation to continually refine a probabilistic estimate of the model count for a formula, so that each XOR-streamlined query yields as much information as possible. We implement this approach, with an approximate probability model, as a wrapper around an off-the-shelf SMT solver or SAT solver. Experimental results show that the implementation is faster than the most similar previous approaches which used simpler refinement strategies. The technique also lets us model count formulas over floating-point constraints, which we demonstrate with an application to a vulnerability in differential privacy mechanisms.

I. INTRODUCTION
Model counting is the task of determining the number of satisfying assignments of a given formula. Model counting for Boolean formulas, #SAT, is a standard model-counting problem, and it is a complete problem for the complexity class #P in the same way that SAT is complete for NP. #P is believed to be a much harder complexity class than NP, and exact #SAT solving is also practically much less scalable than SAT solving. #SAT solving can be implemented as a generalization of the DPLL algorithm [15]. and a number of systems such as Relsat [4], CDP [5], Cachet [32], sharpSAT [35], DSHARP [29] and countAntom [7] have demonstrated various optimization techniques. However, not surprisingly given the problem's theoretical hardness, such systems often perform poorly when formulas are large and/or have complex constraints.
Since many applications do not depend on the model count being exact, it is natural to consider approximation algorithms that can give an estimate of a model count with a probabilistic range and confidence. Some approximate model counters include ApproxCount [38], SampleCount [22], MiniCount [25], ApproxMC [10], ApproxMC-p [37] and ApproxMC2 [11] . In this paper we build on the approximation technique of XOR streamlining [23], which reduces the number of solutions of a formula by adding randomlychosen XOR (parity) constraints. In expectation, adding one constraint reduces the model count by a factor of 2, and k independent constraints reduce the model count by 2 k . If a formula with extra constraints has n > 0 solutions, the original formula likely had about n · 2 k . If the model count after constraints is small, it can be found with a few satisfiability queries, so XOR streamlining reduces approximate model counting to satisfiability. However to have an automated system, we need an approach to choose k values when the model count is not known in advance.
One application of approximate model counting is measuring the amount of information revealed by computer programs. For a deterministic computation, we say that the influence [30] is the base-two log of the number of distinct outputs that can be produced by varying the inputs, a measure of the information flow from inputs to outputs. Influence computation is related to model counting, but formulas arising from software are more naturally expressed as SMT (satisfiability modulo theories) formulas over bit-vectors than as plain CNF, and one wants to count values only of output variables instead of all variables. The theory of arithmetic and other common operations on bounded-size bit-vectors has the same theoretical expressiveness as SAT, since richer operations can be expanded ("bit-blasted") into circuits. But bit-vector SMT is much more convenient for expressing the computations performed by software, and SMT solvers incorporate additional optimizations. We build a system for this generalized version of the problem which takes as input an SMT formula with one bit-vector variable designated as the output, and a specification of the desired precision.
Our algorithm takes a statistical estimation approach. It maintains a probability distribution that reflects an estimate of possible influence values, using a particle filter consisting of weighted samples from the distribution. Intuitively the mean of the distribution corresponds to our tool's best estimate, while the standard deviation becomes smaller as its confidence grows. At each step, we refine this estimate by adding k XOR constraints to the input formula, and then enumerating solutions under those constraints, up to a maximum of c solutions (we call this enumeration process an exhaust-up-to-c query [30]). At a particular step, we choose k and c based on our previous estimate (prior), and then use the query result to update the estimate for the next step (posterior). The update from the query reweights the particle filter points according to a probability model of how many values are excluded by XOR constraints. We use a simple binomial-distribution model which would be exact if each XOR constraint were fully independent. Because this model is not exact, a technique based only on it does not provide probabilistic soundness, even though it performs well practically. So we also give a variant of our technique which does produce a sound bound, at the expense of requiring more queries to meet a given precision goal.
We implement our algorithm in a tool SearchMC 1 that wraps either a bit-vector SMT solver compatible with the SMT-LIB 2 standard or a SAT solver, and report experimental results. SearchMC can be used to count solutions with respect to a subset of the variables in a formula, such as the outputs of a computation, the capability that Klebanov et al. call projected model counting [37], and Val et al. call subset model counting [36]. In our case the variables not counted need not be of bit-vector type. For instance this makes SearchMC to our knowledge the first tool that can be used to count models of constraints over floating-point numbers (counting the floating-point bit patterns individually, as contrasted with computing the measure of a subset of R n as in the work of Chistikov et al. [13]). We demonstrate the use of this capability with an application to a security problem that arises in differential privacy mechanisms because of the limited precision of floating-point values.
Compared to ApproxMC2 [11] and ApproxMC-p [37], concurrently-developed approximate #SAT tools also based on XOR streamlining, our technique gives results more quickly for the same requested confidence levels.
In summary, the key attributes of our approach are as follows: • Our approximate counting approach gives a two-sided bound with user-specified confidence. • Our tool inherits the expressiveness and optimizations of SMT solvers. • Our tool gives a probabilistically sound estimate if requested, or can give a result more quickly if empirical precision is sufficient.

II. BACKGROUND
XOR Streamlining. The main idea of XOR streamlining [23] is to add randomly chosen XOR constraints to a given input formula and feed the augmented formula to a satisfiability solver. One random XOR constraint will reduce the expected number of solutions in half. Consequently, if the formula is still satisfiable after the addition of s XOR constraints, the original formula likely has at least 2 s models. If not, the formula likely has at most 2 s models. Thus we can obtain a lower bound or an upper bound with this approach. There are some crucial parameters to determine the bounds and the probability of the bounds and they need to be carefully chosen in order to obtain good bounds. However, early systems [23] did not provide an algorithm to choose the parameters. Influence. Newsome et al. [30] introduced the terminology of "influence" for a specific application of model counting in quantitative information-flow measurement. This idea can capture the control of input variables over an output variable and distinguish true attacks and false positives in a scenario of malicious input to a network service. The influence of input variables over an output variable is the log 2 of the number of possible output values. Exhaust-up-to-c query. Newsome et al. also introduced the terminology of an "exhaust-up-to-c query", which repeats a satisfiability query up to some number c of solutions, or until there are no satisfying values left. This is a good approach to find a model count if the number of solution is small. Particle Filter. A particle filter [17] is an approach to the statistical estimation of a hidden state from noisy observations, in which a probability distribution over the state is represented non-parametrically by a collection of weighted samples referred to as particles. The weights evolve over time according to observations; they tend to become unbalanced, which is corrected by a resampling process which selects new particles with balanced weights. A particle filtering algorithm with periodic resampling takes the following form: 1. Sample a number of particles from a prior distribution. 2. Evaluate the importance weights for each particle and normalize the weights. 3. Resample particles (with replacement) according to the weights. 4. The posterior distribution represented by the resampled particles becomes the prior distribution to next round and go to step 2.

III. DESIGN
This section describes the approach and algorithms used by SearchMC. It is implemented as a wrapper around an offthe-self bit-vector satisfiability solver that supports the SMT-LIB2 format [3]. It takes as input an SMT-LIB2 formula in a quantifier-free theory that includes bit-vectors (QF BV, or an extension like QF AUFBV or QF FPBV) in which one bit-vector is designated as the output, i.e. the bits over which solutions should be counted. (For ease of comparison with #SAT solvers, SearchMC also has a mode that takes a Boolean formula in CNF, with a list of CNF variables designated as the output.) SearchMC repeatedly queries the SMT solver with variations of the supplied input which add XOR constraints and/or "blocking" constraints that exclude previously-found solutions; based on the results of these queries, it estimates the total number of values of the output bit-vector for which the formula has a satisfying assignment.
SearchMC chooses fruitful queries by keeping a running estimate of possible values of the model count. We model the influence (log 2 of model count) as if it were a continuous quantity, and represent the estimate as a probability distribution over possible influence values. In each iteration we use the current estimate to choose a query, and then update the estimate based on the query's results. (At a given update, the most recent previous distribution is called the prior, and the new updated one is called the posterior.) As the algorithm runs, the confidence in the estimate goes up, and the best estimate changes less from query to query as it converges on the correct result. Each counting query SearchMC makes is parameterized by k, the number of random XOR constraints to add, and c, the maximum number of solutions to count. The result of the query is a number of satisfying assignments between 0 and c inclusive, where a result that stops at c means the real total is at least c. Generally a low result leads to the next estimate being lower than the current one and a high result leads to the estimate increasing. Section III-A describes the process of updating the probability distribution, and then Section III-B gives the details of the algorithms that use it.

A. Updating distribution and confidence interval
We here explain the idea of how we compute a posterior distribution over influence, where both the prior and posterior are represented by particles. Suppose we have a formula f with a known influence log 2 N , and add k XOR random constraints to the formula. If we simulate checking the satisfiability of this augmented formula f for different XOR constraints, we can estimate a probability of sat/unsat on f . We expand this idea by applying exhaust-up-to-c approach to f . We count the number of satisfying assignments n up to c and generate the distributions for each number of satisfying assignments (where n = c means that the number of satisfying assignments is in fact c or more). Thus under an assumption on the true influence of a formula, we can estimate the probabilities of each number of satisfying assignments based on k. By collecting these probabilities across a range of influence, we obtain a probability distribution over influence for an unknown formula assumed to have less than a maximum bits of influence. Under the idealized assumption that each XOR constraint is completely independent, adding k XOR constraints will leave each satisfying assignment alive with probability 1/2 k . For any particular set of n ≥ 0 satisfying assignments remaining out of an original N , the probability that exactly those n solutions will remain is the product of 1/2 k for each n and 1 − (1/2 k ) for each of the other N − n. Summing the total number of such sets with a binomial coefficient, we can approximately model the probability of exactly n solutions remaining as: For the case when the algorithm stops looking when there might still be more solutions, we also want an expression for the probability that the number of solutions is n or more. We compute this straightforwardly as one minus the sum of the probabilities for smaller values: We use XOR constraints that contain each counted bit with probability one half, and are negated with probability one half. (This is the same family of constraints used in other recent systems [10], [37], [11]. Earlier work [23] suggested using constraints over exactly half of the bits, which have the same expected size, but less desirable independence properties.) Our binomial probability model is not precise in general, because these XOR constraints are 3-independent, but not r-independent for r ≥ 4. When N ≥ 4, some patterns among solutions (such as a set of four bitvectors whose XOR is all zeros) lead to correlations in the satisfiability of XOR constraints, and in turn to higher variance in the probability distribution without changing the expectation. This effect is relatively small, but we lack an analytic model of it, so we compensate by slightly increasing the confidence level our tool targets compared to what the user originally requested.
This probability model lets us simulate the probability of various query results as a function of the unknown formula influence. We use this model as a weighting function for each particle and resample particles based on each particle's weight value. Then, we estimate a posterior distribution from sampled particles that have all equal weights. For instance, given a prior distribution over the influence sampled at 0.1 bit intervals, we can compute a sampled posterior distribution by counting and re-normalizing just the probability weights that correspond to a given query result value n. From the estimated posterior distribution, the mean µ and the standard deviation σ are computed. Hence, the µ is our best possible answer as our algorithm iterates and σ shows how much we are close to the true answer. Sequentially, the posterior distribution will be the next round's prior distribution and for use in the very first step of the algorithm we also implement a case of the prior distribution as uniform over influence.
Next we compute a confidence interval (lower bound and upper bound) symmetrically from the mean of the posterior distribution even though the distribution is not likely to be symmetrical. There are several ways to compute the confidence interval but the difference of the results is negligible as the posterior distribution gets narrower. Therefore, we used a simple way to compute the confidence interval: a half interval from the left side of the mean and another half from the right side.

B. Algorithm
We present our main algorithm SearchMC that runs automatically and always gives an answer with a given confidence interval. The pseudocode for algorithm SearchMC is given as Algorithm 1. Our algorithm takes as input a formula f , a desired confidence level CL (0 < CL < 1), a confidence level adjustment α (0 ≤ α < 1), a desired range size thres and an initial prior distribution InitDist. f contains a set of bit-vector variables and bit-vector operators. We can obtain a confidence interval at a confidence level for a given mean and standard deviation. A confidence level CL is a fraction parameter specifying the probability with which the interval should contain the true answer, for example, 0.95 (95%) or 0.99 (99%). As we described above, the soundness of our model has not been proved hence a confidence level sometimes needs to be adjusted by the adjustment value α to meet the actual confidence level requirement. If α = 0, we do not adjust the input confidence level. In our experiments, we used α = 0.5, which empirically leads SearchMC to select large enough ranges. Finding the most appropriate value for α is future work. Our algorithm terminates when the length of our confidence interval is less than or equal to a given non-negative parameter thres. This parameter determines the amount of running time and there is a trade-off. If thres value is small, it gives a narrow confidence interval, but the running time would be longer. If the value is large, it gives a wide confidence interval, but a shorter running time. Our tool can choose any initial prior distribution InitDist represented by particles. For example, if we do not have any knowledge of input formula, it might be better to start with a uniform distribution over 0 to a number of output variables. If we have a small knowledge that the true influence is less than 64, a uniform distribution over 0 to 64 might perform better.
Algorithm 1 SearchMC(f , thres, CL, α, InitDist) The width of the output bit-vector of f 3: prior ← InitDist Initial distribution 4: δ ← width 5: while δ > thres do 6: c, k ← ComputeCandK(prior, width) 7: nSat ← MBoundExhaustUpToC(f, width, k, c) In this case, c is effectively infinite 6: k ← 0 No constraints 7: end if 8: return c, k Variables. There are several variables: prior, post, width, k, c, nSat, U B, LB and δ. prior represents a prior distribution by sampled particles with corresponding weights. In one iteration, we obtain the updated posterior distribution post with resampled particles based on our probabilistic model as described in Section III-A. The posterior becomes the prior distribution for the next iteration. While our algorithm is in the loop, it keeps updating post. width is the width of the output bit-vector of an input formula f , which is an initial upper bound for the influence since the influence cannot be more than the width of the output bit-vector. k is a number of random XOR constraints and c specifies the maximum number of solutions for the exhaust-up-to-c query. We obtain c and k using the ComputeCandK function shown as Algorithm 2 and discussed below. nSat is a number of solutions from the exhaust-up-to-c query. MBoundExhaustUpToC runs until it finds the model count exactly or c solutions from formula f with k random XOR constraints. U B and LB are variables to store an upper bound and a lower bound of the current model count approximation with a given confidence level as we describe in Section III-A. δ is the distance between the upper bound and lower bound. This parameter determines whether our algorithm terminates or not. If δ is less than or equal to our input value thres, our algorithm terminates. If not, it runs again with updated post until δ reaches the desired range size thres. An extreme case k = 0 denotes that our guess is equivalent to the true model count.
In this case, we print out the exact count and terminate the algorithm.
Functions. To motivate the definition of the function ComputeCandK, we view an exhaust-up-to-c query as analogous to measuring influence with a bounded-length "ruler." Suppose that we reduce the expected value of the model count by adding k XOR constraints to f . Then, we can use the "length-(log 2 c) ruler" to measure the influence starting at k and this measurement corresponds to the result of an exhaust-up-to-c query: the length-(log 2 c) ruler has c markings spaced logarithmically as illustrated in Figure 1. Each iteration of the algorithm chooses a location (k) and length (c) for the ruler, and gets a noisy reading on the influence as one mark on the ruler. Over time, we want to converge on the true influence value, but we also wish to lengthen the rule so that the finer marks give more precise readings. Based on this idea, we have the ComputeCandK function to choose the length of and starting point of the ruler from a prior distribution. Then, we run an exhaust-upto-c query and call Update described in Section III-A to update the distribution based on the result of the query.
The pseudocode for algorithm update is described as Algorithm 3. A prior distribution prior and a posterior distribution post are represented as a set of sampled particles (influences). We sampled 500 particles for each update function call. Once we have the updated distribution, we can find out the interval of a given confidence level.
Since we observe that our running σ represents how much we are close to the true answer, we use a rational function to satisfy the condition that c increases as σ decreases (i.e., we get more accurate result as c increases).
The k value denotes where to put the ruler. We want to place the ruler where the expected value of the prior distribution lies near the middle of the ruler hence our expected value is in the range of the ruler with high probability. Therefore, we subtract the half length of the ruler ( 1 2 log 2 c) from the expected value µ and then use the floor function to the x t ← prior t x is a list of sampled particles (influences) 3: if nSat < c then Updating each weight of each particle 4: w t = P r =nSat (2 xt , k) value because k has to be an nonnegative integer value. The expected value always lies in the right-half side of the ruler by using the floor function. However, it is not essential which rounding function is used. Note that there might be a case where k becomes negative. If this happens, we set k = 0 and c = ∞, because our expected value is so small that we can run the solver exhaustively to give the exact model count. The formula for c is motivated by the intuition that the spacing between two marks near the middle of the ruler should be proportional to the standard deviation of the the probability distribution, to ensure that a few different results of the query are possible with relatively high probability; the spacing between the two marks closest to 1 2 . Setting this equal to σ, solving for c, and taking the ceiling gives line 3 of Algorithm 2.

C. Probabilistic sound bounds
We have found that the binomial model performs well for choosing a series of queries, and it yields an estimate of the remaining uncertainty in the tool's results, but because the binomial models differs in an hard-to-quantify way from the true probability distributions, the bounds derived from it do not have any associated formal guarantee. In this section we explain how to use our tool's same query results, together with sound though less-precise bounds, to compute a probabilistically sound lower and upper bound on the true influence. As a trade-off, these bounds are usually not as tight as our tool's primary results.
Lemma 1. Let n = |Σ|, ∈ [0, 1], k ∈ N, k ≤ log 2 (|F | · 2 /(r · 3 √ e)) and h ∈ H(n,k,r) a randomly chosen strongly r-universal hash function. It holds: They called pivot as what we call c from an exhaust-upto-c query and pivot = 2·r·(1+ )· 3 √ e 2 . Since 0 < < 1 and r = 3, c (pivot) should be greater than 17 to make the lemma true with a probability of at least 0.86 (1 − e −3/2 0.86). In SearchMC's iterations, given c and k, we can compute value to estimate the bounds. Therefore, when c is greater than 17 from our tool's iteration, we can compute a lower and upper bound such that the true influence is within the bounds with a probability of at least 0.86.

IV. EXPERIMENTAL RESULTS
In this section, we present our experimental results. All our experiments were performed on a machine with an Intel Core i7 3.40Ghz CPU and 16GB memory. Our main algorithm is implemented with a Perl script and Update function is implemented in a C program called by the main script. Our algorithm can be applied to both SMT formulas and CNF formulas. We have tested a variety of SAT solvers and SMT solvers, and our current implementation specifically supports Cryptominisat2 [33], [34] for CNF formulas and Z3 [16] and MathSAT5 [14] for SMT formulas. We modified Cryp-tominisat2 to count the number of solutions over specified variables. For pure bit-vector SMT formulas, our tool also supports eagerly converting the formula to CNF first and then using CNF mode. (We implement the conversion using the first phase of the STP solver [21], [2] with optimizations disabled and a patch to output the SMT-to-CNF variable mapping.) Performing CNF translation eagerly gives up the benefit of some (e.g., word-level) optimizations performed by SMT solvers, but it can sometimes be profitable because it avoids repeating bit-blasting, and allows the tool to use a specialized multiple-solutions mode of Cryptominisat.

A. Comparison with ApproxMC2 and ApproxMC-p
We run our algorithm with a set of DQMR (Deterministic Quick Medical Reference) benchmarks [1] and ISCAS89 benchmarks [6] converted to CNF files by TG-Pro [12] The DQMR benchmarks are medical diagnosis Bayesian networks and each problem is given by a two layer bipartite network graph for diseases and symptoms. There is an edge between the disease and the symptom if a disease may result a symptom. In these benchmarks, the numbers of diseases and symptoms are varied from 50 to 100 and each symptom   ApproxMC2 and ApproxMC-p are state-of-the-art approximate #SAT solvers which we describe in more detail in Section V. We used Cryptominisat2 as the back-end solver with all the tools for fair comparison. For the parameters for the tools, we set CL = 0.86 (confidence level), α = 0.5 (confidence level adjustment factor) and thres = 1.7 (desired interval length). As described above, SearchMCsound gives correct bounds with a probability of at least 0.86. Since the desired confidence level for ApproxMC2 is 1 − δ, it can achieve a 86% confidence level by setting a parameter δ = 0.14 which corresponds to our parameter CL = 0.86. The length of the interval for ApproxMC2 is computed as log 2 (|f |×(1+ ))−log 2 (|f |×(1/(1+ ))) = 1.7 hence we can obtain the interval length 1.7 by setting a parameter = 0.8, corresponding to our parameter thres = 1.7. Computing the interval for ApproxMC-p is a little different. The length of the interval for ApproxMC-p is log 2 (|f | × (1 + )) − log 2 (|f | × (1 − )) = 1.7 hence we can obtain the interval length 1.7 by setting a parameter = 0.53.
Note that SearchMC increases the c value of an exhaust-upto-c query as it iterates while the corresponding ApproxMC2 and ApproxMC-p parameters are fixed as a function of (72 and 46, respectively) in this experiment. Also, we set an initial prior to be a uniform distribution over 0 to 64 bits for SearchMC. We tested 122 benchmarks (83 DQMRs and 39 ISCAS89s). SearchMC and ApproxMC-p were able to solve 106 benchmarks (83 DQMR and 23 ISCAS89) and ApproxMC2 was able to solve 50 benchmarks (35 DQMRs and 15 ISCAS89s) within 2 hours. The benchmarks that were solved completely by the other tools were also solved completely by SearchMC. Figures 2a and 2d are based on the benchmarks that were solved by all tools. Figure 2a compares the quality of lower bounds and upper bounds computed by SearchMC-sound, ApproxMCp and ApproxMC2. The influence bounds are the computed bounds minus the true influence. Filled markers and empty markers represent reported lower bounds and upper bounds, respectively. Since ApproxMC2 computes the bounds conservatively by using more queries, ApproxMC2's intervals are more closely centered on the true influence, and it outperforms the requested 86% confidence level. By comparison, our tool and ApproxMC-p compute the bounds more aggressively, and in a few cases the true result is just outside the reported interval (visible as empty markers below the dotted line), though this still occurs somewhat less often than the 14% implied by the confidence level. SearchMC-sound tends to give tighter bounds than ApproxMCs since it stops when the interval length becomes less than thres, while the interval lengths for ApproxMCs are fixed by a parameter . Figure 2b shows another perspective on the trade-off between performance and error. We selected a single benchmark and varied the parameter settings of each algorithm, measuring the absolute difference between the returned answer and the known exact result. We include results from running ApproxMC2 with parameter settings outside the range of its soundness proofs (shown as "disallowed" in the plot), since these settings are still empirically useful, and SearchMC makes no such distinction. From this perspective the tools are complementary depending on one's desired performance-error trade-off. The results from all the tools improve similarly with configurations that use more queries, but SearchMC performs best at getting more precise results from a small number of queries. About 100 satisfiability queries is a minimum for ApproxMC2 to give any results. On the other hand our current SearchMC and SearchMCsound implementation's results improve as using more queries, the error varies over different runs compared to the stable error of ApproxMC2.
We also compare the running-time performance with Ap-proxMCs and show the running-time performance comparison on our 122 benchmarks in Figure 2c. Since ApproxMCp refined the formulas of ApproxMC, it used a smaller number of queries than ApproxMC2. SearchMC can solve all the benchmarks faster than ApproxMCs with 86% confidence level. SearchMC-sound performs faster than ApproxMC-p even SearchMC-sound computes its confi-dence interval similarly to ApproxMC-p. The SearchMC's and SearchMC-sounds average running times are 41.79 and 84.89 seconds, compared to an average of 127.8 for ApproxMC-p. ApproxMC2 requires an average of 281.25 seconds just for the subset of benchmarks it can complete.
We also compare the number of SAT queries on the benchmarks for all the tools in Figure 2d. The average number of SAT queries for SearchMC, SearchMCsound ApproxMC-p and ApproxMC2 is about 29.11, 79.75, 1257.19 and 4613.72 queries, respectively. Again this average for ApproxMC2 is based on 50 benchmarks and others are based on 106 benchmarks.

B. Floating Point / Differential Privacy Case Study
As an example of model counting with floating point constraints, we measure the security of a mechanism for differential privacy which can be undermined by unexpected floating-point behavior. The Laplace mechanism achieves differential privacy [20] by adding exponentially-distributed noise to a statistic to obscure its exact value. For instance, suppose we wish to release a statistic counting the number of patients in a population with a rare disease, without releasing information that confirms any single patient's status. In the worst case, an adversary might know the disease status of all patients other than the victim; for instance the attacker might know that the true count is either 10 or 11. If we add random noise from a Laplace distribution to the statistic before releasing it, we can leave the adversary relatively unsure about whether the true count was 10 or 11, while preserving the utility of an approximate result. A naive implementation of such a simple differentially private mechanism using standard floating-point techniques can be insecure because of a problem pointed out by Mironov [28]. For instance if we generate noise by dividing a random number in [1, 2 31 ] by 2 31 and taking the logarithm, the relative probability of particular floating point results will be quantized compared to the ideal probability, and many values will not be possible at all. If a particular floating point number could have been generated as 10 + noise but not as 11 + noise in our scenario, its release completely compromises the victim's privacy.
To measure this danger using model counting, we translated the standard approach for generating Laplacian noise, including an implementation of the natural logarithm, into SMT-LIB 2 floating point and bit-vector constraints. (We followed the log function originally by SunSoft taken from the musl C library, which uses integer operations to reduce the argument to [ √ 2/2, √ 2), followed by a polynomial approximation.) A typical implementation might use double-precision floats with an 11-bit exponent and 53-bit fraction, and 32 bits of randomness, which we abbreviate "53e11, 2 32 ", but we tried a range of increasing sizes. We measured the total number of distinct values taken by 10 + noise as well as the size of the intersection of this set with the 11 + noise set.
The results and running time are shown in Table I. We ran SearchMC with a confidence level of 80%, a confidence level adjustment of 0.5 and a threshold of 1.0; the SMT solver was MathSAT 5.3.13 with settings recommended for floating-point constraints by its authors. We use one random bit to choose the sign of the noise, and the rest to choose its magnitude. The sign is irrelevant when the magnitude is 0, so the expected influence for n bits of randomness is log 2 (2 n − 1). SearchMC's 80% confidence interval included the correct result in all cases. The size of the intersections is small enough that SearchMC usually reports an exact result, and always much less than the total set of noise values, confirming that using this algorithm and parameter setting for privacy protection would be illadvised. The running time increases steeply as the problem size increases, which matches the conventional wisdom that reasoning about floating-point constraints is challenging. But because floating-point SMT solving is a young area, there is significant scope for future solvers to improve the technique's performance.

V. RELATED WORK
In this section, we summarize related work on model counting. We categorize model counting techniques into three areas: exact model counting, randomized approximate model counting and non-randomized approximate model counting. We also introduce some applications of model counting techniques for security and privacy purposes.
Exact model counters give the exact number of solutions for a formula but don't perform well as the size of a problem increases. Approximate model counters have been proposed to resolve the scalability challenges. When a formula has a large number of solutions, it is often sufficient to provide rough estimates instead of the exact model counts, especially when this can be done much faster. Randomized approximate model counting techniques are likely to use random sampling and a SAT solver to produce a probabilistic result, which can often be argued to provide a lower bound and/or upper bound with high probability. In contrast, non-randomized approximate model counting techniques generally use approximations that are more efficient but do not provide probabilistic bounds.

A. Exact model counting
Some of the earliest Boolean model counters used the DPLL algorithm [15] for counting the exact number of solutions. Birnbaum et al. [5] formalized this idea and introduced an algorithm for counting models of propositional formulas. Based on this idea, Relsat [4], Cachet [32] sharpSAT [38] and DSHARP [29] showed improvements by using several optimizations. Relsat uses component analysis in which the model count of a formula is the product of the model count of each sub-formula (component). Cachet shows optimizations by combining component caching and clause learning. sharpSAT introduces an improved caching technique to reduce the space requirement compared to Cachet. DSHARP [29] is a CNF → d-DNNF compiler for efficient reasoning and uses sharpSAT as a back-end. The major contribution of countAntom [7] is techniques for parallelization, but it provides state-of-the-art performance even in single-threaded mode.
Phan et al. [31] encode a full binary search for feasible outputs in a bounded model checker. This approach is precise, but requires more than one call to the underlying solver for each feasible output. Specifically, it recursively calls the solver by adding a bit constraint for finding a single satisfying assignments as DPLL-based search. This is most similar to the exhaust-up-to-c approach that requires n + 1 queries to find n feasible outputs. This search tree approach is useful when the program verification system does not expose the underlying logical representation or when the used solver cannot generate models. Klebanov et al.. [24] perform exact model counting for quantitative informationflow measurement, with an approach that converts C code to a CNF formula with bounded model checking and then uses exact #SAT solving. They explore both exhaustive enumeration and the existing DSHARP and sharpSAT tools, but only counting distinct values of the output variables. Val et al. [36] integrate a symbolic execution tool more closely with a SAT solver by using techniques from SAT solving to prune the symbolic execution search space, and then perform exact model counting restricted to an output variable.
However, precisely counting solutions is a #P-complete problem and these exact model counters typically only work well with small-sized problems or ones with only simple constraints. For many practical problems, it is infeasible to count the exact number of solutions in a reasonable amount of time.
B. Randomized approximate model counting Randomized approximate model counting techniques perform well on many kinds of a formula for which finding single solutions is efficient. Also, there can sometimes be a smooth trade-off chosen between computational effort and the precision of results. However, this solving is still relatively expensive, so research to get the best precision for a given cost is still important. Wei and Selman [38] introduced ApproxCount which uses near-uniform sampling to estimate the true model count but it can significantly over-estimate or underestimate if the sampling is biased. SampleCount [22] improves this sampling idea and gives a lower bound with high probability by using a heuristic sampler. MiniCount [25] is based on a framework to compute an upper bound under statistical assumptions which is that counting the number d of branching decisions (except unit propagations and failed branches) can be used to estimate the total number of solutions by setting a variable to true or false randomly. Specifically, they showed that the expected value of d is not lower than log 2 of the true model count. By estimating the expected value, they can obtain the upper bound of the true model count. Also, they observed that d often has a distribution that is close to a normal distribution. Thus, the expected value can be easily computed under this assumption rather than computing low and high values of d. This only guarantees the upper bound and they used a different method to compute a lower bound.
MBound [23] is an approximate model counting tool that gives probabilistic bounds on the model counts by adding randomly-chosen parity constraints as XOR streamlining. Chakraborty et al. [10] introduced ApproxMC, an approximate model counter for CNF formulas, which automated the choice of XOR streamlining parameters. The ApproxMC algorithm, in our terminology, starts by fixing c and a total number of iterations based on the desired precision and confidence of the results. In each iteration ApproxMC searches for an appropriate k value, adds k XOR random constraints, and then performs an exhaust-up-to-c query on the streamlined formula and multiplies the result by 2 k . It stores all the individual estimates as a multiset and computes its final estimate of the model count as the median of the values. The original ApproxMC sequentially increases k in each iteration until it finds an appropriate k value. An improved algorithm ApproxMC2 [11] uses galloping binary search and saves a starting k value between iterations to make the selection of k more efficient. Other recent systems that build on ApproxMC include SMTApproxMC [9] and ApproxMC-p [37]. SMTApproxMC proposes word-level constraints based on modular arithmetic instead of bit-level XOR constraints, however these will not likely provide com-parable performance until SMT solvers implement modulararithmetic Gaussian elimination. ApproxMC-p implements projection (counting over only a subset of variables), and also gives more efficient formulas for parameter selection.
ApproxMC2, developed concurrently, is the system most similar to SearchMC: its binary search for k plays a similar role to our converging µ value. However SearchMC also updates the c parameter over the course of the search, leading to fewer total queries. ApproxMC, ApproxMC2, and related systems choose the parameters of the search at the outset, and make each iteration either fully independent (ApproxMC) or dependent in a very simple way (ApproxMC2) on previous ones. These choices make it easier to prove the tool's probabilistic results are sound, but they require a conservative choice of parameters. By comparison SearchMC's approach of maintaining a probabilistic estimate explicitly at runtime means that its iterations are not at all independent: instead our approach is to extract the maximum guidance for future iterations from previous ones, to allow the search to converge more aggressively.
The runtime performance of SearchMC, like that of ApproxMC (2), is highly dependent on the performance of SAT solvers on CNF-XOR formulas. Some roots of the difficulty of this problem have been investigated by Dudek et al. [18], [19].

C. Non-randomized approximate model counting
Non-randomized approximate model counting using techniques similar to static program analysis is generally faster than randomized approximate model counting techniques, and such systems can give good approximations for some problem classes. However, they cannot provide a precision guarantee for arbitrary problems, and it is not possible to give more effort to have more refined results.
Castro et al. [8] compute an upper bound on the number of bits about an input that are revealed by an error report. They measure the entropy loss of an error report by computing the number of bits revealed by subsets of path conditions first and then combining these partial results to get the final result. Meng and Smith [27] use two-bit-pattern SMT entailment queries to calculate a propositional overapproximation and count its instances with a model counter from the computer algebra system Mathematica. Luu et al. [26] propose a model counting technique over an expressive string constraint language. Their tool computes the bounds on the cardinality of the valid string set and uses generating functions for reasoning about the cardinality of string sets.

D. Applications: Security and Privacy
Various applications of model counting have been proposed for security and privacy purposes. For instance Castro et al. [8] use model counting and symbolic execution approaches to measure leaking private information from bug reports. They compute an upper bound on the amount of private information leaked by a bug report and allow users to decide on whether to submit the report or not. Newsome et al. [30] show how an untrusted input affect a program and introduce a family of technique for measuring influence which can be applicable to x86 binaries. They applied the XOR streamlining and the exhaust-up-to-c techniques separately, but did not combine the two techniques or describe an algorithm for choosing XOR streamlining parameters.

VI. FUTURE WORK
Closing the gap between the performance of SearchMC and SearchMC-sound is one natural direction for future research. On one hand, we would like to explore techniques for asserting sound probabilistic bounds which can take advantage of the results of all of SearchMC's queries. At the same time, we would like to find a model of the number of solutions remaining after XOR streamlining that is more accurate than our current binomial model, which should improve the performance of SearchMC. Another future direction made possible by the particle filter implementation is to explore different prior distributions, including unbounded ones. For instance, using a negative exponential distribution over influence as a prior would avoid the any need to estimate a maximum influence in advance, while still starting the search process with low-k queries which are typically faster to solve.

VII. CONCLUSION
In sum, we have presented a new model counting approach SearchMC using XOR streamlining for SMT formulas with bit-vectors and other theories. We demonstrate our algorithm that adaptively maintains a probabilistic model count estimate based on the results of queries. Our tool computes a lower bound and an upper bound with a requested confidence level, and yields results more quickly than previous systems.


Fully Dynamic Algorithm for Top-$k$ Densest Subgraphs

Given a large graph, the densest-subgraph problem asks to find a subgraph with maximum average degree. When considering the top-$k$ version of this problem, a na\"ive solution is to iteratively find the densest subgraph and remove it in each iteration. However, such a solution is impractical due to high processing cost. The problem is further complicated when dealing with dynamic graphs, since adding or removing an edge requires re-running the algorithm. In this paper, we study the top-$k$ densest-subgraph problem in the sliding-window model and propose an efficient fully-dynamic algorithm. The input of our algorithm consists of an edge stream, and the goal is to find the node-disjoint subgraphs that maximize the sum of their densities. In contrast to existing state-of-the-art solutions that require iterating over the entire graph upon any update, our algorithm profits from the observation that updates only affect a limited region of the graph. Therefore, the top-$k$ densest subgraphs are maintained by only applying local updates. We provide a theoretical analysis of the proposed algorithm and show empirically that the algorithm often generates denser subgraphs than state-of-the-art competitors. Experiments show an improvement in efficiency of up to five orders of magnitude compared to state-of-the-art solutions.

INTRODUCTION
Finding a subgraph with maximal density in a given graph is a fundamental graph-mining problem, known as the densest-subgraph problem. Density is commonly de ned as the ratio between number of edges and vertices, while many other de nitions of density have been used in the literature [7,28,30,31]. e densest-subgraph problem has many applications, for example, in community detetion [11,14], event detection [2], link-spam detection [18], and distance query indexing [1].
In applications, we are o en interested not only in one densest subgraph, but in the top-k. e top-k densest subgraphs can be vertex-disjoint, edge-disjoint, or overlapping [6,17]. Di erent objective functions and constraints give rise to di erent problem formulations [6,17,32]. In this work, we choose to maximize the sum of the densities of the k subgraphs in the solution. In addition, Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for pro t or commercial advantage and that copies bear this notice and the full citation on the rst page. Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permi ed. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior speci c permission and/or a fee. Request permissions from permissions@acm.org.  Figure 1b shows the top-3 densest subgraphs a er the arrival. e objective is to design an algorithm that can e ciently maintain the densest subgraphs while keeping the number of updates very low, in this case updating only the vertices in red.
we seek a solution with disjoint subgraphs.
is version of the problem is known to be NP-hard [6].
To complicate the ma er, most real-world graphs are dynamic and rapidly changing. For instance, Facebook users are continuously creating new connections and removing old ones, thus changing the network structure. Twi er users produce posts at a high rate, which makes old posts less relevant. Given the dynamic nature of many graphs, here we focus on a sliding-window model which gives more importance to recent events [4,12,13]. Finding the topk densest subgraphs in a sliding window is of interest to several real-time applications, e.g., community tracking [33], event detection [26], story identi cation [2], fraud detection [8], and more. We assume the input to the system arrives as an edge stream, and seek to extract the k vertex-disjoint subgraphs that maximize the sum of densities [6].
A naïve solution involves executing a static algorithm for the densest-subgraph problem k times, while removing the densest subgraph in each iteration. However, such a solution is impractical as it requires to execute the algorithm k times for each update. An alternative solution to our problem is to use a dynamic densest subgraph algorithm in a pipeline manner, where the output of an algorithm instance serves as input to the following one. In this case, the graph and the instances of the algorithm are replicated independently across k instances of the algorithm, resulting in a high memory and processing cost.
In this paper, we propose a fully-dynamic algorithm that nds an approximate solution. e proposed algorithm follows a greedy approach and updates the densities of the subgraphs connected to vertices a ected by edge operations (addition and removal). e algorithm is e ciently designed based on key properties of densest subgraphs, and it is competitive against other recent dynamic algorithms [9,15,24].
First, our algorithm relies on the observation that only highdegree vertices are relevant for the solution. As many natural graphs have a heavy-tailed degree distribution, the number of highdegree vertices in a graph is relatively smaller than the number of low-degree ones. is simple observation enables pruning a major portion of the input stream on-the-y. Second, the vertices that are part of a densest subgraph are connected strongly to each other and weakly to other parts of the graph. is enables independently maintaining and locally updating multiple subgraphs. Figure 1 provides an example which demonstrates this intuition. e algorithm tracks multiple subgraphs on-the-y with the help of a newly de ned data structure called snowball. ese subgraphs are stored in a bag, from which the k subgraphs with maximum densities are extracted. e algorithm runs in-place, and does not require multiple copies of the graph, thus making it memorye cient. e one-pass nature of the algorithm allows extracting top-k densest subgraphs for larger values of k.
We provide a theoretical analysis of the proposed algorithm, and show that the algorithm guarantees 2-approximation for the rst densest subgraph (k = 1) while providing a high-quality heuristic for k > 1 compared to other solutions. Experimental evaluation shows that our algorithm o en generates denser subgraphs compared to the state-of-the-art algorithms, due to the fact that it maintains disconnected subgraphs separately. In addition, the algorithm provides improvement in runtime up to three to ve orders of magnitude compared to the state-of-the-art. In summary, we make the following contributions: • We study the top-k densest vertex-disjoint subgraphs problem in the sliding-window model. • We provide a brief survey on adapting several algorithms for densest subgraph problem for the top-k case. • We propose a scalable fully-dynamic algorithm for the problem, and provide a detailed analysis of it. • e algorithm is open source and available online, together with the implementations of all the baselines. 1 • We report a comprehensive empirical evaluation of the algorithm in which it signi cantly outperforms previous state-ofthe-art solutions by several orders of magnitude, while producing comparable or be er quality solutions.

PRELIMINARIES
In this section, we present our notation, revisit basic de nitions, and formulate the top-k densest subgraphs problem. Consider an undirected graph G = (V , E) with n = |V | vertices and m = |E| edges. e neighborhood of ∈ V is de ned as N ( ) = {u | ( , u) ∈ E}, and its degree as d( ) = |N ( )|. For a subset S ⊆ V we de ne E(S) to be the set of edges whose both endpoints are in S, and G(S) = (S, E(S)) the subgraph induced by S. e internal degree of a vertex with respect to S ⊆ V is de ned by d S ( ) = |N ( ) ∩ S |. 1 h ps://github.com/anisnasir/TopKDensestSubgraph Finally, for a subset of vertices S ⊆ V we de ne its density ρ S by Note that the density of any subgraph is equal to half of its average internal degree.
De nition 2.1 (Densest subgraph). Given an undirected graph G = (V , E), the densest subgraph S * is a set of vertices that maximizes the density function, i.e., We say that an algorithm A computes an α-approximation of the densest subgraph if A computes a subset S ⊆ V such that ρ S ≥ 1 α ρ S * , where S * ⊆ V is the densest subgraph of G. Next we introduce other concepts related to densest subgraph: graph core, core decomposition, and induced core subgraph of a vertex.
De nition 2.2 (j-core). Given an undirected graph G = (V , E) and an integer j, a j-core of G is a subset of vertices C ⊆ V so that each vertex ∈ C has internal degree d C ( ) ≥ j, and C is maximal with respect to this property.
where each C i is a j-core for some j.
De nition 2.4 (Core number). Given a core decomposition V = C 0 ⊇ C 1 ⊇ . . . ⊇ C ⊇ ∅ of a graph G = (V , E), the core number κ( ) of a vertex is the largest j such that ∈ C and C is a j-core.
By overwriting notation, the core number κ(C) of a core C is the largest j for which C is a j-core.
Additionally, we use κ S ( ) to denote the core number of a vertex in the subgraph induced by S. e largest core (or main core) of a subgraph of G(S) = (S, E(S)) is denoted by C (S), while the main core of G is simply denoted by C .
Note that the density of a j-core is at least j /2, as each vertex in the core has degree at least j and each edge is counted twice. is observation implies that the main core of a graph is a 2-approximation of its densest subgraph. L 2.5. Consider the core decomposition of a graph G, i.e., C 1 ⊆ C 2 ⊆ . . . ⊆ C . e maximum core C is 2-approximation of the densest subgraph of G [21].

P
. Let S * be the densest subgraph of G having density ρ S * . Every vertex in G(S * ) has degree at least ρ S * ; otherwise a vertex with degree smaller than ρ S * can be removed to obtain an even denser subgraph.
us, S * is a ρ S * -core. Given the core decomposition of G, we know that ρ C ≥ 1 2 κ(C ). We want to show that ρ C ≥ 1 2 ρ S * . Assume otherwise, i.e., ρ C < 1 2 ρ S * . en κ(C ) < ρ S * . It follows that S * is a higher-order core than C , a contradiction. e concept of a core subgraph I( ) induced by a vertex [23,27] is also pertinent to our analysis.
De nition 2.6 (Induced core subgraph). Given a graph G = (V , E) and a vertex ∈ V , the induced core subgraph of vertex , denoted by I( ), is a maximal connected subgraph containing the vertex s.t. the core number of all the vertices in I( ) is equal to κ( ).
In other words, the induced subgraph contains all vertices that are reachable from and have the same core number κ( ).
All previous de nitions apply to static graphs. Let us now focus on dynamic graphs. In particular, we consider processing a graph in the sliding window edge-stream model [13]. According to this model, the input to our problem is a stream of edges. e edge e i is the i-th element of the stream. Equivalently, we say that edge e i has timestamp i. A sliding window W t (x), de ned at time t and of size x, is the set of all edges that arrive between e t −x +1 and e t , For each edge e i = (u, ), we consider that u and appear at time i, and we use V t (x) to denote the set of vertices that appear in a length-x sliding window at time t. e graph in a length-x sliding window at time t is then de ned to be G t (x) = (V t (x),W t (x)).
We are now ready to formally de ne the problem that we consider in this paper, i.e., nding the top-k densest subgraphs in sliding window. We rst de ne the problem in a static se ing.
De nition 2.7. Given an undirected graph G = (V , E) and an integer k > 0, the top-k densest subgraphs of G is a set of k disjoint maximal set of vertices S = {S 1 , . . . , S k } that maximize the sum of its densities: As already shown by Balalau et al. [6], the problem de ned aboce is NP-hard, for any k > 1. e problem we consider in this paper is the following. P 2.8. Given a graph stream {e i } and a sliding window length x, maintain the top-k densest subgraphs ρ k (S) of the graph G t (x), at any given time t.

BACKGROUND
In this section we present a brief review over several algorithms for nding dense subgraphs. Additionally, we discuss how these methods can be used for solving Problem 2.8. Densest subgraph in static graphs. Finding the densest subgraph according to the density de nition (1) can be solved in polynomial time. An elegant solution involving reduction to the minimumcut problem was given by Goldberg [20]. As the fastest algorithm to solve the minimum-cut problem runs in O(nm) time [25], Goldberg's algorithm is not scalable to large graphs.
is algorithm iteratively removes the vertex with the lowest degree in each iteration, until le with an empty graph. Among all subgraphs considered during this vertex-removal process, the algorithm returns the densest. e time complexity of this greedy algorithm is O(m + n).
Bahmani et al. [5] propose a MapReduce version of the greedy algorithm, with approximation ratio 2(1 + ϵ), while making O(log 1+ϵ n) passes over the input graph. Core decomposition in static graphs. e core decomposition of a graph G is the process of identifying all cores of G, as de ned in 2.3. Batagelj and Zaversnik [7] propose a linear-time algorithm to obtain the core decomposition. e algorithm rst considers the whole graph and then repeatedly removes the vertex with the smallest degree. e core number κ( ) of a vertex is set equal to the degree of at the moment that is removed from the graph. Densest subgraph in evolving graphs. ere is a growing body of literature on nding dense subgraphs in evolving graphs [9,15,19,24]. We focus mainly on the deterministic algorithm for densest subgraph in evolving graphs. For instance, Epasto et al. [15] propose an e cient algorithm for computing the densest subgraph in the dynamic graph model [16]. eir work assumes that edges are inserted into the graph adversarially but deleted randomly. Even though the algorithm can, in practice, handle arbitrary edge deletions, its approximation guarantees hold only under the random edge-deletion assumption. e algorithm is similar to the one by Bahmani et al. [5], and it provides a 2(1 + ϵ) 6 -approximation of the densest subgraph, while requiring polylogarithmic amortized cost per update with high probability. Core maintenance in evolving graphs. Saríyüce et al. [27] propose the traversal algorithm, for e cient core maintenance. is algorithm identi es a small set of vertices that are a ected by edge updates and processes these vertices in linear time in order to maintain a valid core decomposition. Li et al. [23] propose an e cient three-stage algorithm for core maintenance in large dynamic graphs. e algorithm maintains a core decomposition of an evolving graph by applying updates to very few vertices in the graph. Once these few vertices have been identi ed, the algorithm computes the correct core numbers via a quadratic operation. Finding top-k densest subgraphs. e problem of nding top-k densest subgraphs has been mainly studied for nding overlapping subgraphs in static graphs [6,17].
Our rst observation is that a set of k dense subgraphs can be obtained from any algorithm that nds the densest subgraph by k repeated invocations. e time complexity of computing a set of k dense subgraphs in this manner is simply the runningtime complexity of the densest-subgraph algorithm multiplied by k. From a practical point of view, all the static algorithms mentioned are not in-place algorithms, and thus require copying the whole graph for processing. Furthermore, when a vertex or edge is added or deleted from the graph, the whole k dense subgraph computation has to be repeated. e second observation is that, by using the algorithm of Epasto et al. [15], we can obtain a set of k dense subgraphs by running k instances of the fully-dynamic algorithm in a pipeline manner. e idea is to run k instances of the algorithm in which the output of each instance i ∈ {0 . . . k − 1} is fed into the next (i + 1) instance as a removal operation. e pipeline version of the algorithm requires keeping k copies of the input graph and an additional O(kn) size space for bookkeeping. Note that the output of each instance of the pipeline might cascade, which requires updating the vertices in all the instances. In particular, vertices that cease to be part of solution in upstream instances need to be added in downstream instances. Likewise the vertices that become part of densest subgraphs in upstream instances need to be removed from a downstream instances. e modi cation of the algorithm, as discussed above, is expensive in terms of memory, as it requires replicating the graph and the algorithm's structures k times. In addition, running and maintaining k parallel instances makes the algorithm compute-intensive.
Finally, to maintain top-k densest subgraphs in evolving graphs, we can leverage algorithms for core decomposition maintenance [23,27]. by leveraging Lemma 2.5, us, the idea is to nd and maintain top-k disjoint maximum cores. In order to maintain such cores we run a single instance of the algorithm by Saríyüce et al. [27] or Li et al. [23] that maintains the core number of all the vertices in the graph. We then extract the top-k densest subgraphs by: (i) extracting the main core, (ii) removing the vertices in the main core and updating the core number for rest of the vertices, and (iii) repeating the steps until k subgraphs are extracted.

ALGORITHM
e main idea of our algorithm is to maintain and update multiple dense subgraphs online. ese subgraphs are candidates for the top-k densest subgraphs. However, maintaining multiple subgraphs for fully-dynamic streams requires answering two interesting questions: (i) how to reduce the search space of the solution, and (ii) how to split the whole graph into subgraphs.
To answer the aforementioned questions, we make two observations. First, since dense subgraphs are formed by relatively highdegree vertices, one can nd dense subgraphs by keeping track of these high-degree vertices only. Second, these subgraphs can be updated locally upon edge updates, without a ecting the other parts of the graph.
Based on these observations, we develop an algorithm that reduces the solution space by considering only high-degree vertices, and divides the whole graph into smaller subgraphs, each representing a dense subgraph. e top-k densest subgraphs among the candidate subgraphs provide a solution for Problem 2.8.
We begin by designing an algorithm to nd the densest subgraph (top-1) and then we extend it to nd the top-k densest subgraphs. Our algorithm might not be the most e cient solution for the (top-1) densest-subgraph problem per se, but it provides e cient outcomes when extended to solve the top-k densest-subgraph problem.
We start by de ning some properties of the densest subgraph that we leverage in our algorithm.
is lemma holds according to the de nition of optimal density. In an optimal solution, each vertex has degree larger than or equal to ρ S * . Otherwise, removing the vertex from the subgraph will increase the average degree, and thus the density, of the subgraph. Given Lemma 4.1, at any time t, the densest subgraph S * t of graph G t contains only vertices that have degree d( ) ≥ d S * t ( ) ≥ ρ S * t . en, given G t and ρ S * t , we want to compute the densest subgraph a er the addition of a new vertex u V t at time t + 1.
Let d(u) be the degree of vertex u ∈ V t +1 and S * t +1 be the densest subgraph at time t + 1. For simplicity, assume that the graph G t +1 is connected. According to Lemma 4.1, for any vertex u to be the part of the densest subgraph, its internal degree satis es d S * t +1 (u) ≥ ρ S * t +1 . As vertex u is added to the graph the new density is always greater, i.e., ρ S * t +1 ≥ ρ S * t . erefore, for vertex u to be the part of densest subgraph, the degree of vertex u should satisfy d(u) ≥ ρ S * t . erefore, if the degree of vertex u is lower than the ρ S * t , it cannot be part of the densest subgraph ρ S * t +1 and can be ignored. Now, considering the case when d(u) ≥ ρ S * t . Adding the vertex u to densest subgraph will update the density: We also know that ρ S * t +1 ≥ ρ S * t , which means that Using these properties, we ignore the vertices of the new edge that have degree lower than the current estimate of the density.
Further, we are interested in nding the main core in the remaining subgraph of high-degree vertices, as it represents a 2approximation of the densest subgraph according to Lemma 2.5. To this end, we propose a new data structure that relies on Lemma 4.1, the snowball.

Snowball
A snowball D is an incremental data structure that stores a strongly connected subgraph, which maintains the following invariants: • e core number κ D ( ) of each vertex ∈ D inside a snowball is equal to the main core (C (D)) of the snowball.
• All the vertices in the snowball are connected. ese invariants ensure that all the vertices in the snowball have the same core number, which is the main core of the snowball by de nition. A snowball maintains these invariants while handling the following graph update operations: a) adding/removing a vertex, and b) adding/removing an edge.

Bag of snowballs
e high-degree vertices in the graph are assigned to a snowball. As these vertices are not strongly connected, they might end up in di erent snowballs. We store each of these disconnected snowballs in a data structure called the bag, denoted by B.
e bag ensures that each snowball is vertex disjoint. Further, the bag provides an additional operation: extracting the densest snowball among the set of snowballs. e density of the extracted snowball is the maximal density, which is the threshold separating the high-degree vertices from the low-degree ones. We denote this estimate of the maximal densityρ S * . e bag is a supergraph which contains a set of snowballs and all the edges between the snowballs. We maintain all the core  Figure 2: Example showing that the bag requires maintaining the core number of the vertices. Initially, the bag contains two snowballs with core number 2, i.e., D 1 and D 2 . Consider the arrival of the edges shown in red. e greedy assignment of the edges might skip creating a new snowball with core number 3, using the four nodes in the middle.
numbers of the nodes in the bag by leveraging a core decomposition algorithm (see Section 3). e core number of each node in the bag is used to ensure that each node has the maximum possible core number. Figure 2 provides an example explaining one of the issues that may arise. In the example, the bag contains two snowballs, however, it is possible to produce a new snowball with a larger core number. Next, we de ne the algorithms to update this data structure upon graph updates.

Addition operations
Vertex addition: As discussed in Section 2, the updates appear in the form of an edge stream. Here, we de ne the vertex addition algorithm that acts as a helper for edge addition. e algorithm is triggered when at least one of the endpoints of a new edge is a high-degree vertex. In particular, there are two cases to consider: 1) the bag already contains the high-degree vertex, and 2) the bag does not contain the high-degree vertex. In both cases, the goal is to add the new vertex to one of the snowballs (if needed).
Algorithm 1 de nes the algorithm for vertex addition. For the rst case, the algorithm scans the bag to nd the snowball that contains the vertex and returns it. For the second case, the algorithm rst identi es the candidate snowballs, then it assigns the vertex to one of the candidate snowballs. e candidate snowballs are the ones having the main core number smaller than or equal to the internal degree of the new vertex (κ u (D) ≥ C (D)). Among the candidate snowballs, the new node is assigned to the snowball with maximum internal degree d u (D), breaking ties randomly.
Once a vertex is added to a snowball, the core number of the snowball may increase. is change requires removing the vertices with core number lower than the main core of the snowball. is procedure can be implemented e ciently in linear time by sorting the vertices based on their degree similar to bin sort. 2 Veri cation. Due to the greedy assignment of vertex to the snowball, it is possible that the vertex ends up not having the highest possible core number. For example, Figure 3 shows an example where the greedy assignment does not result in optimal solution. erefore, a er addition, the algorithm ensures that the core number of the snowball, where the new vertex is added, equals the core number of the new vertex in the graph. e algorithm veri es that the core number of the added vertex by comparing it with the core number of the vertex in the bag. Note that the bag represents the supergraph containing all the snowballs and edges between the snowballs. If the core number within the bag is  S * ← ∅ 3: r epeat ← f alse 4: if (d (u) ≥ρ S * ) then 8: B (u) 9: r epeat ← t rue 10: while r epeat larger than the one in the snowball, the algorithm merges all the snowballs in the induced subgraph of newly added vertex. As all the vertices in the induced subgraph have the same core number, merging them ensures creating a larger snowball. 3 We leverage the core decomposition algorithm by Saríyüce et al. [27] for the implementation.

Algorithm 3 Fix Main Core
Given the bag B, the algorithm ensures that B contains the main core of the graph within one of the snowballs a er the vertex addition.

P
. Let us assume that at time t the bag contains the main core of the graph. Now, we need to show that at time t + 1, a er the node addition, the bag still contains the main core of the graph. In general, vertex addition method is called whenever there is an edge addition. e only way for the new vertex to a ect the main 3

P
. e proof for all the cases, other than the case when both the vertices of the new edge are assigned to two di erent snowballs, is similar to the vertex addition algorithm. erefore, we consider the case when both end vertices of the added edge are added to two di erent snowballs. For this case, we rely on the graph in the bag. We check if both vertices are in the same core graph in the bag, and x the core number of the vertices if they belong to the same induced subgraph. is ensures creating the graph with the highest core number.

Removal operations
Vertex removal: Similarly to the addition operations, we rst dene the procedure for removing a vertex from the bag of snowballs. e vertex removal method is used as a subroutine for the edge removal process. A vertex is only removed from the bag when its degree becomes lower than the maximal density. erefore, according to Lemma 2.5, the removed vertex cannot be part of the main if ((d (u) <ρ S * ) and (d ( ) <ρ S * )) then 3: return 4: if ((d (u) <ρ S * ) and (d (u) + 1 ≥ρ S * )) then 5: V (u) 6: return 7: if ((d ( ) <ρ S * ) and (d ( ) + 1 ≥ρ S * )) then 8: V (v) 9: return 10: for D i ∈ B do 11: if (D i ∩ (u, ) ∅) then 12: D i ← D i \(u, ) 13: core. e algorithm removes the vertex from the snowball within a bag without doing any other operation.
Edge removal: Now we turn our a ention to edge deletion, which follows the same pa ern as edge addition. Again, we leverage the bag to ensure that there exist a snowball with a core number equal to the main core of the graph. Algorithm 5 shows the algorithm for edge deletion. e bag does not require any update when either one of the vertices is low-degree or if both the vertices belong to two di erent snowballs. erefore, we consider the case when one of the vertices lie at the boundary of high-degree vertices. at is, the edge deletion moves the vertex from the high-degree to low-degree. In this case, the algorithm only requires removing the vertex from the bag, without performing any other operations. e interesting case is where both vertices of the deleted edge are high-degree and belong to the same snowball. In this case, the algorithm removes the edge from the snowball. Further, it veri es and updates (if needed) the core number of the vertices a ected by the update in the snowball. Lastly, edge deletion might reduce the maximal density, and thus require adding to the bag new vertices whose degree is now greater than the new maximal density. T 4.4. Given the bag containing a snowball with the same core number as the main core of the bag, a er the edge deletion, Algorithm 5 maintains the main core of the graph in one of the snowballs inside the bag.

P
. An edge removal a ects the bag of snowballs only when both the vertices corresponding to the removed edge belong to the same snowball. In this case, edge removal might reduce the core number of all the vertices in the snowball. e algorithm ensures that the vertices in the snowball have their maximum possible core number by comparing their core number in the snowball with the core number in the bag. erefore, by verifying and xing the core numbers, the algorithm maintains the main core of the graph in one of the snowballs inside the bag.

Fully-dynamic top-k densest subgraphs
Now that we have a fully dynamic algorithm for nding the densest subgraph in sliding windows, we move our a ention to the top-k densest-subgraph problem.
Letρ S * ≥ ρ 1 ≥ . . . ≥ ρ z−1 represent the densities of the z subgraphs in the bag, where ρ S * is the density of the densest subgraph. e bag contains the vertices that have a degree greater than the densityρ S * . To ensure that the bag contains at least k subgraphs, we modify the algorithm to keep all the vertices with degree greater than ρ k −1 . e only modi cation required is to replace ρ S * by ρ k −1 in Algorithm 2, Algorithm 4 and Algorithm 5. Further, we leverage the priority queue to extract ρ k−1 from the bag of snowballs. is simple modi cation ensures that the bag contains at least k subgraphs and enables accessing the top-k densest subgraphs in the sliding window model. Note that the new de nition of high-degree vertices is related to the density of the k-th top densest subgraph. e modi ed algorithm guarantees that the bag contains the vertices with a degree greater than the ρ k −1 a er any graph update. ese graph updates include both edge additions and removals. It is necessary for the algorithm to consider edge additions, as they might a ect the value of ρ k −1 by merging multiple subgraphs. Similarly, edge deletions have to be considered, as they might a ect the value of ρ k −1 by reducing the density of any of the top-k densest subgraphs, merging multiple subgraphs, or spli ing a subgraph.
As the algorithm ensures keeping the main core in the bag, it guarantees 2-approximation for the rst densest subgraph (k = 1) while providing a high-quality heuristic for k > 1 compared to other solutions. We provide an example in Figure 4 for the top-k densest subgraph algorithm by expanding on the Figure 1a. We conclude this section with the theorem that provides a bound for our proposed algorithm. ese bounds can be generalized for any algorithm that adapts to a solution for densest subgraph problem for the top-k case (see section 3 for examples). Given an integer k, the bag contains a set of subgraphs that provides 2k-approximation of the top-k densest-subgraph problem.

P
. We know that the graph does not contain any subgraph with density greater than the optimal density (ρ S * ). erefore, we know that the sum of densities for the top-k densest-subgraph problem is upper bounded by k × ρ S * .
Further, the bag containsρ S * , which provides a 2-approximation of the densest subgraph.
is implies that the sum of densities of top-k densest subgraph in the bag ≥ 1 2 ρ S * . Pu ing above two observations together, we can clearly see that the bag provides a solution that is a 2k-approximation for the problem.

Data structure
e proposed algorithm requires accessing the neighborhood information of every node. Speci cally, we are interested in performing three queries on a given vertex: a) extract its degree, b) extract all its neighbors, and c) given densityρ S * , extract all the vertices with degrees greater than the densityρ S * . Vertex map: To answer the rst two queries, we need to store the neighborhood information for all vertices. We store the information in a hashmap with keys being the vertex identi ers and values being the neighbors of each vertex. Vertex map allows performing search and update operations in amortized constant time. Degree table: For the third query, we need to order the vertices by their degrees. We use bin sort to order the vertices by their degrees, which enables extracting the vertices with degrees greater than the densityρ S * in constant time.

DISCUSSION
Most of the solutions that leverage a static algorithms [5,7,10] require iterating over the entire graph k times upon any update. Comparatively, our algorithm mostly touches a limited region in the graph for any updates, which makes our algorithm perform signi cantly faster. e top-k densest subgraph algorithm that adapts to Epasto et al. [15] requires replicating the graph across the k instances. In addition, updates across k instances might cascade and require running the algorithm k times, similarly to the static top-k solutions. Our algorithm outperforms the pipeline version of the top-k algorithm both in terms of memory and computation.
Finally, incremental algorithms for k-core maintenance requires removing top-k densest subgraph upon every update in the graph. Each removal of a densest subgraph requires updating the core number of the remaining vertices. Our algorithm e ciently avoids this removal step by maintaining the subgraphs during execution.
Performance. e proposed algorithm leverages the skew in realworld graphs and ignores a major portion of the input stream on the y. In addition, high-degree vertices are stored as small subgraphs so that each update is o en applied locally on these subgraphs. is design allows our algorithm to operate on just small portions of the graph for each update, rather than iterating on the whole graph. Finally, we use the k-core algorithm [23,27], which allows each high-degree node to update their core number with a complexity independent of the graph size. ese improvements enable our algorithm to perform signi cantly be er than the other state-ofthe-art streaming algorithms. Memory. Our algorithm requires only O(n 2 polylogn) memory, compared to O(kn 2 polylogn) for the top-k version of [15] and O(n 2 polylogn) for the other algorithms discussed in Section 3. Tight Bounds. We show via an example that any algorithm for densest subgraph problem can only produce a k-approximation for top-k densest subgraph problem. Consider a graph G that contains n ·m vertices, for any n ≥ k. e n nodes connect to form a circle. In the circle, there is one edge connecting two non-adjacent vertices. Additionally, each of the n nodes connects to exactly m neighbors. e inner circle of the graph G is the densest subgraph. Removing the densest subgraph leaves the rest of the vertices completely disconnected. Hence, the sum of density will be equal to one. Now, consider the case when the densest subgraph is not removed rst.
In this case, each node in the circle along with its m neighbors create a subgraph of density almost equal to one, for higher values of m. erefore, one can return k such subgraphs with the sum of densities equal to k, which is k times be er than the previous case.

EVALUATION
We conduct an extensive empirical evaluation of the proposed algorithm, and provide comparisons with the existing solutions. In particular, we answer the following questions: Q1: What is the impact on the quality of subgraphs? Q2: What are the gains in performance? Q3: How does the algorithm perform in terms of di erent input parameters?

Experimental setup
Datasets. Table 1 shows the datasets used in the experiments. e datasets are selected due to their public availability. We evaluate all the algorithms in the sliding window model. e number of edges in the sliding window is an input parameter. Metrics. We evaluate the quality and e ciency of the algorithms. We assess the quality by the objective function, i.e., the sum of densities of the subgraphs produced by an algorithm. We evaluate the e ciency of an algorithm by reporting the average update time and the memory usage. e average update time is the average time it takes to move the sliding window.
is includes adding the new edge, removing the oldest edge, and updating the top-k densest subgraphs. We report the memory usage as the average percentage of occupied memory. Algorithms. Table 2 shows the notations used for di erent algorithms. e algorithms by Bahmani et al. [5] and Epasto et al. [15] require an additional epsilon parameter for execution, which provides a trade-o between quality and execution time. Here we use the defaults proposed by the authors. Stream ordering. We consider two commonly used stream ordering schemes [29]: • BFS: e ordering is a result of a breadth-rst search starting from a random vertex. • DFS: e ordering is a result of a depth-rst search starting from a random vertex. Experimental environment. We conduct our experiments on a machine with 2 Intel Xeon Processors E5-2698 and 128GiB of In this experiment, we compare the quality of the results produced by the di erent algorithm as measured by the objective function. In order to be able to run most of the algorithms, we use the smaller datasets, i.e., AM, DB 1 , YT, and DB 2 . As RL and T R produce same results in terms of quality, we only report the results for one of them. Due to space constraints, we show the results only with the DFS ordering. However, we achieve similar results also with the BFS one. We set the size of the sliding window x = 100k. For the static algorithms, we execute them in micro-batches in which the top-k densest subgraph is recomputed a er 1k edge updates, which gives them a substantial advantage. For BB ϵ and AE ϵ , we use ϵ = 0.01. We set the maximum execution time to 7 days, due to which the AE ϵ is not able to nish on DB 2 and YT. Finally, RL and T R contain an update phase due to the iterative extraction of the top-k densest subgraph. is phase is expensive, as the main core consists of high-degree vertices. To alleviate this issue, we execute the extraction phase every 10k update operations. Figure 5 shows the quality results. Our algorithm GR achieves competitive quality, o en generating denser subgraphs compared to all the other algorithms. For example, for the AM dataset, GR produces subgraphs that are 1.5 times be er than the best algorithm among the state-of-the-art solutions. All the other algorithms produce consistent results across all the datasets. For instance, V B, RL, and T R produce top-k dense subgraphs of lowest quality compared to the other algorithms. is result is caused by the removal of the main core (backbone) of the graph in each of the k iterations. In addition, it shows how the problem considered in this paper, while related, is di erent from a simple k-core decomposition. e results also validate that BB ϵ and AE ϵ provide weaker guarantees on the quality compared to CH . Q2: In this experiment, we turn our a ention to evaluate the efciency of the proposed algorithm, as measured by the average update time and memory usage. We again select several datasets, i.e., AM, DB 1 , YT, DB 2 , and LJ, and execute CH , V B, BB , RL, T R, AE , and GR. We exclude the largest datasets as only a few algorithms are able to handle them. e algorithms are executed with both BFS and DFS ordering of the edge stream. e sliding window size is set to x = 100k edges, and k = 10 unless otherwise speci ed. 6 h ps://github.com/anisnasir/TopKDensestSubgraph   Figure 6 shows the e ciency results (note the logarithmic scale). e proposed algorithm, GR, outperforms all other ones in terms of update time for all the datasets and both ordering schemes. RL and T R are the slowest algorithms. In particular, for DFS, our algorithm achieves a performance gain of ve orders of magnitude. is result is noteworthy as even though our algorithm is dependent on core decomposition, it is still able to beat a naïve application of those algorithms by a wide margin. is di erence is due to the fact that both algorithms require maintaining the core number for all the vertices. Additionally, the extraction of top-k densest subgraph further hampers their e ciency.
Algorithms CH , V B, and BB 0.01 perform unfavorably due to their static nature. In this case, our algorithm is able to achieve more than three orders of magnitude improvement in e ciency.
Algorithm AE 0.01 is best performing among the baselines, and even outperforms GR in one dataset for one speci c ordering (AM with BFS). However, GR still outperforms AE 0.01 by nearly three orders of magnitude in most other cases.
Finally, we observe the overall memory consumption of all the algorithms (reported in Table 3). e memory requirement of our algorithm lies between the static and the dynamic algorithms, i.e., AE 0.01 requires the largest amount of memory, while the static algorithms require the least. ese results are in line with our expectations from the discussion in Section 3. Q3: In this experiment, we evaluate the scalability of our algorithm. First, we execute GR with di erent values of k. Alongside, we report the average update time for BB ϵ algorithm. We choose BB ϵ , rather than AE ϵ , as it can execute in mirco-batches, i.e., top-k densest subgraph every 100 edges. We set the sliding window size x = 1M, and use the YT and DB 2 dataset with DFS ordering. Figure 7 reports the average update time for both algorithms. e average update time of GR remains consistent even for higher values of k, whereas the execution time of BB 0.01 increases with the parameter.
Second, we execute GR with di erent sizes of sliding window, i.e., x = 10k, 100k, 1M, and 10M. We set k = 10 for this experiment. Again, we use the YT and DB 2 datasets with DFS ordering.    Table 4 reports the average update time for di erent con gurations. Increasing the size of the sliding window does not a ect the average update time signi cantly. is result validates our claim that most of the updates are local to the subgraphs, and do not require iterating through the whole graph to extract the top-k densest subgraphs. Finally, we study the performance of our algorithm in on the largest datasets. We select the OT and FR datasets with DFS ordering, and execute the algorithm with a sliding window of x = 100k, with k = 10. Figure 8 shows the result of the experiment. e plot is generated by taking the moving average of the update time and the sum of densities. e average update time of the algorithm mostly remains constant throughout the execution, and our algorithm provides steady e ciency over the stream.
is behavior remains consistent even when the densities are uctuating, as in the case for the OT dataset.

RELATED WORK
Valari et al. [32] were the rst one to study the top-k densest subgraph problem for a stream consisting of a dynamic collection of graphs. ey proposed both an exact and an approximation algorithm for top-k densest subgraph discovery. Similar to our algorithm, the proposed algorithm relies on the core decomposition to provide the approximation guarantees. e top-k densest subgraphs produced by the algorithm are edge-disjoint. Balalau et al. [6] studied the problem of nding the top-k densest subgraph with limited overlap. ey de ned the top-k densest subgraphs as a set of k subgraphs that maximizes the sum of densities, while satisfying an upper bound on the pairwise Jaccard coe cient between the set of vertices of the subgraphs. e problem of nding the top-k densest subgraph as sum of densities was shown to be NP-hard [6] and e cient heuristic was proposed to solve the problem. Further, Galbrun et al. [17] studied the problem of nding the top-k overlapping densest subgraphs and provided constant-factor approximation guarantees.

CONCLUSION
We studied the top-k densest subgraphs problem for graph streams, and proposed an e cient one-pass fully-dynamic algorithm. In contrast to the existing state-of-the-art solutions that require iterating over the entire graph upon update, our algorithm maintains the solution in one-pass. Additionally, the memory requirement of the algorithm is independent of k. e algorithm is designed by leveraging the observation that graph updates only a ect a limited region. erefore, the top-k densest subgraphs are maintained by simply applying local updates to small subgraphs, rather than the complete graph. We provided a theoretical analysis of the proposed algorithm and showed empirically that the algorithm o en generates denser subgraphs than the state-of-the-art solutions. Further, we observed an improvement in performance of up to ve orders of magnitude when compared to the baselines.
is work gives rise to further interesting research questions: Is it necessary to leverage k-core decomposition algorithm as a backbone? Is it possible to achieve stronger bounds on the threshold for high-degree vertices? Can we design an algorithm with a space bound on the size of the bag? Is it possible to achieve stronger approximation guarantees for the problem? We believe that solving these questions will further enhance the proposed algorithm, making it a useful tool for numerous practical applications.