1 Introduction↩︎

Score-based Diffusion Models have presented the deep learning community with a very, if not the most, attractive way of acquiring highly local information about a complex high-dimensional probability distribution [9]. These models work by training a neural network to approximate the score (gradient of its log function) of an evolving probability density function, which, starting out as the distribution to be learned, is gradually turned into Gaussian noise through a diffusion process. Using the learned score function, score-based diffusion models can then reverse the diffusion process by solving a backward diffusion stochastic differential equation (SDE), and, in doing so, sample from [[1]][10][11] or evaluate [[1]][3] the underlying distribution at the user’s request.
Given their success, one may naturally wonder whether cardinality estimation, a task crucial for the ideal performance of query optimizers, can also benefit from diffusion models. Indeed, several learned cardinality estimation models, such as Naru[4] and DQM-D[12], are already performing cardinality estimation by approximating an underlying probability density function at individual points, then numerically integrating it across the queried region using various Monte Carlo schemes[[4]][13]. However, simply grabbing a popular score-based diffusion model off the shelf probably won’t work due to the three reasons below.
First, while conventional models used for image generation excel at approximating distributions with thousands of dimensions, their sampling process often take up to tens of seconds even with the help of advanced TPUs[14] and time-saving techniques like consistency models[15] or DPM Solver[16], which is often longer than the time it costs to actually execute the query.
Second, existing works on diffusion models have never touched on the problem of integrating the probability density function across a given region, a task useless for image generation, yet crucial for cardinality estimation.
Most importantly, image generation and cardinality estimation require different magnitudes of precision: while making it five times more likely to generate a specific image gives a generative model its unique art style, predicting five thousand tuples to exist in a region where only one thousand actually does gives a query optimizer an optimization disaster.
In this paper, we aim to overcome the aforementioned limitations by presenting Accelerated Diffusion Cardest(ADC), a “downsized” diffusion model much more suited for cardinality estimation, structured according to the framework in Figure321.

Here, our major contributions are:

1. Downsizing Score Matching. We find that for distributions with moderate dimensionality, the inherent structure and properties of diffusion processes can be used to construct much more lightweight score matching networks that achieve the same level of precision. That is, score matching is best done using the data prediction model when \(t\) is large, and a modified noise prediction model when \(t\) is small. We modify traditional noise prediction models by constructing a neural network with three output modules, normalized using the functions \(\frac{1}{\sigma^2(t)}, \frac{1}{\sigma(t)}\), and \(1\), respectively. This modification is proven to significantly enhance precision and, to the best of our knowledge, has not been proposed in previous works.

2. Simplifying Likelihood Calculation. Existing works studying diffusion models [3] provide a formula for calculating the pointwise log-density of a probability distribution by integrating its score function across time and space. We convert this formula into a fast, well grounded algorithm which combines Quasi Monte Carlo (for integration across space) with adaptive stepsize Midpoint Rule (for integration across time), enabling accurate estimations in less than 1 ms per point.

3. GMM for Cardinality Estimation. In order to estimate the selectivity of ranged queries, we turn to the importance sampling algorithm \(\int_V p(x)dx=\int_Vq(x)dx\mathrm{E}_{q(x|x\in V)}\frac{p(x)}{q(x)}\), where the problem boils down to constructing a suitable predictor \(q\). Here, we propose training a Gaussian Mixture Model (GMM) for \(q\), as conditional sampling, pointwise density calculation, and integration across hyperrectangles all prove very easy for GMMs. Numerical experiments find that GMMs not only cooperate very well with diffusion models, but also prove to be efficient standalone estimators. This is especially true for queries with large selectivity or great volumes, an advantage we harnessed in ADC+ by training a decision tree to identify queries "easy" enough for the GMM \(q(x)\) to handle independently.

4. Experimental Evaluation Highlights. We test ADC and ADC+ against five state-of-the-art learned models (DeepDB[7], Naru[4], LW-NN, LW-Tree[8], MSCN[6]), strictly adhering to well-established benchmarks for dataset choice and workload design[2]. We find that on real-world datasets, ADC+ is comparable to Naru and almost always beats MSCN, DeepDB, LW-Tree, and LW-NN, using less than 360KB of storage space and with a latency at most half that of Naru’s. Further experiments using synthetic datasets show that ADC works best on relations where strong synergistic interactions exist between attributes, quantifiable by a large and positive interaction information (i.e. one can, and can only, significantly narrow down the range of one attribute by knowing several others). In this case, no other learned model performs as well as ADC+, and Naru even loses to ADC+ on 95%th and 99%th Q-error by a more than 10 \(\times\) margin.

2 Background, Notations, and Related Work↩︎

2.1 Cardinality Estimation↩︎

The goal of cardinality estimation is to rapidly and accurately estimate the number of rows returned by a given query without actually executing it. This essay, in particular, studies a more specific and highly influential subcase presented by [17] and [8] (implied in their problem formulation and experiment dataset choice, though LW-Tree and LW-NN versions supporting categorical attributes are later developed [2]), in which the underlying database is a single relation \(\mathrm{R}\) with continuous attributes \(\{A_1,...,A_n\}\), and the query \(\mathrm{Q}\) is of the form
SELECT count(*) FROM \(\mathrm{R}\)
WHERE \(a_{k_1}<A_{k_1}<b_{k_1}\)
AND \(a_{k_2}<A_{k_2}<b_{k_2}\)
AND ...
AND \(a_{k_i}<A_{k_i}<b_{k_i}\),
i.e. the queried region is always a hyper-rectangle.
Current approaches to cardinality estimation can be very roughly categorized into learned and non-learned methods [2]. Traditional non-learned methods, such as small-scale sampling, multidimensional histograms [18], Bayesnet [19], and kernel density estimators [17], are usually more time-efficient and relatively robust to constant updates [2], but much less accurate than learned models. Therefore, they are of more use when the data distribution is simple, when the database is frequently updated, or when the precision requirements are not strict.
Learned models, on the other hand, usually conduct cardinality estimation in one of the following two ways [2]. Regression models, such as MSCN [6] and LW-tree [8], learn from already-conducted queries by conveying them into feature vectors, then constructing a regression model to match these vectors with the true cardinality of their corresponding queries. Joint distribution models, meanwhile, assume that data points are distributed according to an unknown probability density function, and thus turn cardinality estimation into a numerical integration problem: namely, calculating the density function at individual points, then summing them up using a predefined integration scheme. For example, estimators like Naru [4] and DQM-D [12] calculate the joint distribution function at each point by factorizing it into conditional distributions \[P(A_1,...,A_n)=\prod_{i=1}^dP(A_i|A_1,...,A_{i-1}),\] then sum them up using the Monte Carlo scheme \[\int_V p(x)dx=\int_V g(x)dx\cdot\mathrm{E}_{g(x|x\in V)}\frac{p(x)}{g(x)},\] with the predictor, \(g(x)\), constructed in multiple stages using sequential sampling (Naru), VEGAS [13] (DQM-D), or other sampling algorithms.
Given that estimating and sampling from complex, high-dimensional distributions is exactly what score-based diffusion models do best, as is proven by their groundbreaking success in image, text, and video generation [20], it is no surprise that in the field of joint distribution estimators, diffusion models have the full potential to walk under the spotlight.

2.2 Score-Based Diffusion Models↩︎

Score-based diffusion models are the generative models that made AI drawing and video generation a reality, thanks to pioneering work by Song, Dickstein and Kingma [1], Song and Dhariwal [15], Ho, Jain, and Abbeel [14], Leobacher and Pillichshammer [11], and many others. For the sake of convenience, all mathematical notations related to diffusion models shall be documented in Appendix A. The meaning of each notation will also be explained upon its first appearance.
As mentioned before, given a set of data points satisfying an unknown distribution \(p_0\) in the sample space \(\mathbb{R}^d\), score-based diffusion models attempt to learn and sample from \(p_0\) in two steps: the Forward Process and the Backward Process[1].
The Forward Process gradually injects white noise into the sample points to create Gaussian kernels around them, slowly turning \(p_0\) into a smooth Gaussian distribution using the SDE \[dx=-\alpha(t)xdt+\beta(t)dw,t\in [0,T],\] where \(dw\) refers to the standard Brownian motion.
Two most common choices for \(\alpha\) and \(\beta\) are \(\alpha(t)=0,\beta (t)=1\) (the VE scheme) and \(\alpha(t)=1,\beta(t)=\sqrt{2}\) (the VP scheme, which is used in our model)[1]. However, works including [21] and [22] show that we can, in theory, choose any positive function for \(\alpha\) and \(\beta\) without fundamentally changing the model. In fact, for each choice of \(\alpha\) and \(\beta\), it would always hold that ³ \[\begin{align} p_{t,\alpha,\beta}&=(\tau_{k(t)}\circ p_0)*\varphi_{\sigma^2(t)}\\ \nabla \log p_{t,\alpha,\beta}(x)&=k(t)\nabla \log p_{\sigma^2(t)k^2(t),0,1}(k(t)x), \end{align}\] where \(k(t)=e^{\int_0^t\alpha(s)ds}\), \(\sigma^2(t)\) satisfies \(\frac{d\sigma^2}{dt}=\beta(t)-2\alpha(t)\sigma^2(t)\) with initial condition \(\sigma^2(0)=0\), \(\tau\) in \(\mathbb{R}^d\) is the scaling operator \((\tau_k\circ f)(x)=k^df(kx)\), and \(\varphi_{\sigma^2(t)}\) in \(\mathbb{R}^d\) is the probability density function of the normal distribution \(\mathcal{N}(0,\sigma^2(t)I_d)\).
The Backward Process, in turn, generates a sample point by first sampling a starting point \(\tilde{x}(T)\sim\mathcal{N}(0,\sigma^2(T)I_d)\approx p_{T,\alpha,\beta}(x)\), then use it as an initial value to solve the reverse diffusion SDE [1] \[dx=-[\alpha(t)x+\beta^2(t)\nabla \log p_t(x)]dt+\beta(t)d\hat{w}\] (where \(\hat{w}\) denotes a standard Wiener process in the reverse time direction), or its associated probability flow ordinary differential equation (ODE) \[\frac{dx}{dt}=-[\alpha(t)x+\frac{1}{2}\beta^2(t)\nabla \log p_t(x)].\] Here, \(\nabla \log p_t(x)\), the only unknown term in both equations, is called the score function. To approximate it, diffusion models train a neural network \(s_\theta(x,t)\) to minimize the score prediction error \[\begin{align} \int_0^T\mathrm{E}_{p_{t,appr}(x)}[\beta^2(t)\left\|\nabla \log p_{t,appr}(x)-s_\theta (x,t)\right\|_2^2]dt, \end{align}\] in a process called score matching, with \[\begin{align} p_{t,appr}(x)=\frac{1}{N}\sum_{i=1}^N\mathcal{N}(\frac{x_i}{k(t)},\sigma^2(t)I_d) \end{align}\] being our "best guess" for \(p_t(x)\) from \(N\) i.i.d. samples.
The accuracy of score-based diffusion models has been verified both experimentally and theoretically. In fact, the KL divergence between \(p_0(x)\), the original underlying distribution, and \(p_{0,\theta(x)}\), which \(\tilde{x}(0)\) is distributed accordingly, would satisfy [3] \[\begin{align} D_{KL}\left(p_0\left(x\right),p_{0,\theta}\left(x\right)\right)\leq D_{KL}\left(p_{T,\alpha,\beta},\mathcal{N}\left(0,\sigma^2(T)I_d\right)\right) \notag \\ +\frac{1}{2}\int_0^T\mathrm{E}_{p_{t,\alpha,\beta}(x)}[\beta^2(t)\left\|\nabla \log p_{t,\alpha,\beta}(x)-s_\theta (x,t)\right\|_2^2]dt, \end{align}\] where the equal sign holds whenever \(s_\theta (x,t)\) is the exact score function of a diffusion process (not necessarily \(\nabla \log p_{t,\alpha,\beta}(x)\)).
Noting that the above formula did not require how "similar" \(p_0\) must be from \(p_{0,\theta}\), replacing \(p_0(x)\) with \(\mathcal{N}(x_0,\epsilon^2I_d)\), a normal distribution centered at an arbitrary point \(x_0\), and letting \(\epsilon \to 0\) produces a more interesting formula [3] \[\log p_{0,\theta}(x_0)\approx\mathrm{E}_{p_{0T}(x|x_0)}\log \varphi_{\sigma^2(T)}(x)+\int_0^T n\alpha(t)dt+\frac{1}{2}\int_0^T\beta^2(t)\] \[\mathrm{E}_{p_{0t(x|x_0)}}[\left\|\nabla \log p_{0t}(x|x_0)-s_\theta(x,t)\right\|_2^2-\left\|\nabla \log p_{0t}(x|x_0)\right\|_2^2]dt,\] where \(p_{0t}(x|x_0)=\mathcal{N}(\frac{x_0}{k(t)},\sigma^2(T))\) is the distribution derived by perturbing \(\delta_{x_0}\) for time \(t\) under the same perturbation scheme. Again, the equal sign would exactly hold when \(s_\theta (x,t)\) is the score function of some diffusion process, a condition that would almost be met if \(s_\theta\) can match \(\nabla \log p_{t,\alpha,\beta}(x)\) sufficiently well.
The above formula lays the cornerstone for our model, allowing us to estimate a probability density function at any point, and hence conduct cardinality estimation, by training \(s_\theta(x,t)\) and integrating it across time and sample space. In section 3 and section 4, we shall design the score estimator and density estimator that can finish the task in a matter of milliseconds.

3 The Score Estimator↩︎

3.1 Model Choice↩︎

The score estimator trains a neural network, \(s_\theta(x,t)\), to approximate the score function \(\{\nabla\log p_t(x)|t\in(0,T)\}\), thus providing the density estimator with the data it needs to calculate pointwise densities. Compared to training a network that directly outputs \(\nabla\log p_t(x)\), works on score-based diffusion models usually find one of the following two methods much more efficient.
These approaches, named the data prediction model [[23]][24] and the noise prediction model [[1]][23] respectively, are both inspired by a corollary of formula (1) that \(x(t)\sim p_t(x)\) can be expressed as the sum of two independent variables \(x_{sig}\) and \(x_{noise}\), where \(k(t)x_{sig}\sim p_0(x)\) and \(\sigma(t)^{-1}x_{noise}\sim\mathcal{N}(0,I_d)\). Therefore, \(\nabla\log p_t(x)\) would equal \[\frac{\int_{\mathbb{R}^d}[\tau_{k(t)}\circ p_0](x_{sig})\nabla\varphi_{\sigma^2(t)}(x-x_{sig})dx_{sig}}{\int_{\mathbb{R}^d}[\tau_{k(t)}\circ p_0](x_{sig})\varphi_{\sigma^2(t)}(x-x_{sig})dx_{sig}}=\mathrm{E}_{p_s(x_{sig}|x)}\frac{x_{sig}-x}{\sigma^2(t)},\] where \[p_s(x_{sig}|x)=[\tau_{k(t)}\circ p_0](x_{sig})\varphi_{\sigma^2(t)}(x-x_{sig})\] denotes the conditional distribution of \(x_{sig}\) for a fixed \(x\).
As a result, the data prediction model works by having \(v_\theta(x,t)\) learn to match \(\mathrm{E}_{p_s(x_{sig}|x)}k(t)x_{sig}\), the normalized conditional expectation of \(x_{sig}\) for a given \(x\), and then output \[\nabla\log p_t(x)\approx s_\theta(x,t)=\frac{v_\theta(x,t)}{\sigma^2(t)k(t)}-\frac{x}{\sigma^2(t)}.\] Similarly, defining \[p_n(x_{noise}|x)=[\tau_{k(t)}\circ p_0](x-x_{noise})\varphi_{\sigma^2(t)}(x_{noise}),\] the noise prediction model has \(v_\theta(x,t)\) learn to match \[\mathrm{E}_{p_n(x_{noise}|x)}\frac{x_{noise}}{\sigma(t)},\] and then output \(\nabla \log p_t(x)\approx s_\theta(x,t)=\frac{v_\theta(x,t)}{\sigma(t)}\).
Conventional works on diffusion models usually stick to one of the two designs throughout the score matching process. However, we find that identifying the best-fit scenario of each model can prove remarkably easy. Indeed, for any choice of \((x,t)\), absolute error from approximating \(\mathrm{E}_{p_s(x_{sig}|x)}x_{sig}\) and \(\mathrm{E}_{p_n(x_{sig}|x)}x_{noise}\) equally affect the error term \(\left\|\nabla \log p_{t,approx}(x)-s_\theta (x,t)\right\|_2^2\), as \(x_{sig}\) and \(x_{noise}\) determine each other by a simple plus-and-minus relation. Thus, predicting the smaller of \(x_{sig}\) and \(x_{noise}\), that is, \(x_{noise}\) when \(t\) is small and \(x_{sig}\) when \(t\) is large, naturally leads to a larger relative error tolerance.
With this in mind, we have ADC’s score estimator break score matching into two parts: approximating \(\nabla \log p_t(x)\) with a data prediction model \(s_{\theta_2,tail}(x,t)\) when \(t\) is above a threshold \(T_{trunc}\), and a noise prediction model \(s_{\theta_1,head}(x,t)\) when \(t\) is below it. To the best of our knowledge, this is the first time such a practice has been adpoted in the field of score-based diffusion models.
Many studies [15] [9] [20] show that \(p_{t,appr}(x)\), when constructed from i.i.d. samples, grows drastically larger and more volatile as \(t\to 0\) and eventually explodes to infinity. As a result, \(s_{\theta_1,head}\) must be somewhat heavy to produce moderately accurate estimations, but \(s_{\theta_2,tail}\) can be very lightweight yet just as precise, especially since when \(p_t\) gets sufficiently close to \(\mathcal{N}(0,\sigma^2(t))\), the data prediction model’s added term, \[-\frac{x}{\sigma^2(t)}=\nabla\log \varphi_{\sigma^2(t)}(x),\] has already completed the bulk of score matching pretty decently. Therefore, our approach of implementing two separate models not only improves precision, but can also speed up computation by calculating the score function with a heavy network only where it is really needed.
Finally, we would like to point out that the extra storage space cost incurred by training two separate models is quite minimal: In all our test scenarios, the extra network used by the data prediction model, \(s_{\theta_2,tail}\), only takes up slightly more than a third of \(s_{\theta_1,head}\)’s storage space.

3.2 Model Architecture Improvement: Introducing QuadNet↩︎

Due to difficulties encountered in training, we modify the noise prediction model based on the observation that the term \(x_{noise}\) might or might not cancel itself out when calculating \(\mathrm{E}_{p_n(x_{noise}|x)}\frac{x_{noise}}{\sigma(t)}\), making the growth pattern of \(\nabla \log p_t(x)\) vary significantly depending on the smoothness of \(p_0(x)\) instead of being fixed to a uniform \(\mathcal{O}(\frac{1}{\sigma(t)})\), as is assumed by the traditional model design.
The theoretical foundation of our proposed modifications is summed up in the theorem below.

The above findings inspire us to present a network structure called Quadnet. That is, for a \(d\)-dimensional dataset, we have our network \(s_{\theta_1,head}\) generate a vector \((v_{\sigma^2},v_\sigma,v_1)\) of length \(3d\), and output \[\nabla\log p_t(x)\approx \frac{v_{\sigma^2}}{\sigma^2(t)}+\frac{v_\sigma}{\sigma(t)}+v_1,\] so that each module effectively captures the behavior of \(\nabla \log p_t(x)\) near one type of region, while correcting the residual error term near the more "explosive" regions.
To illustrate the advantage of Quadnet, we define a toy dataset as of Figure322, and have a Quadnet of shape \([3,10,25,25,10,6]\) compete against three compare group networks of shape \([3,10,25,25,10,6,2]\) in performing the score matching task for small \(t\). Compared to its competitors that possess only the quadratic (\(\frac{1}{\sigma^2(t)}\)), linear (\(\frac{1}{\sigma(t)}\), i.e. a traditional noise prediction model), or constant (i.e. directly approximates \(\nabla \log p_t(x)\)) scaling module, Quadnet converges faster, settles to a lower error, and becomes less prone to falling into saddle points than all of its competitors.

3.3 Mitigating Score Blowup↩︎

Despite the enhancements above, it is widely observed that the score function would still grow to a point beyond the approximation ability of any neural networks when \(t\) becomes too close to zero. Score-based diffusion models used for image generation mitigate this problem by implementing early stopping [15] [20]: that is, training the score matching model (i.e., calculating the score approximation loss) and running the reverse diffusion SDE solver on the interval \((\epsilon,T]\) rather than \((0,T]\), and treating \(x(\epsilon)\) rather than \(x(0)\) as the final sample.
The score estimator in ADC is trained using the same early stopping approach, with \(\epsilon\), the early stopping time, selected among values \(\{\frac{1}{640},\frac{1}{320},\frac{1}{160}\}\) using a trial-and-error approach (starting from \(\frac{1}{640}\) and increasing \(\epsilon\) if the MSE loss proves too big). We notice that the best choice for \(\epsilon\) increases with dimensionality and volume of the data distribution’s typical set, but more work is still needed to determine the best choice for \(\epsilon\) ahead of model construction.
Another natural consequence of score blowup, theoretically coined in the data prediction model case by [9], is that the score approximation error when \(t\) is small always accounts for the bulk of the entire error term. In fact, in our model where \(T\) is set to 3 and \(T_{trunc}\) is set to \(1/8\), \[\int_{\epsilon}^{\frac{1}{32}}\mathrm{E}_{p_t(x)}\left\|\nabla\log p_t(x)-s_\theta(x,t)\right\|_2^2dt\] often contributes to more than 60% of (3)’s size. Therefore, a score estimator’s quality can almost entirely be judged by whether it can perform score matching with enough accuracy on the lower end of the evaluation interval.
With this in mind, ADC speeds up score matching by iterating between exploration epoches and adjustment epoches in its training process. In exploration epoches, the model is told to solely focus on reducing approximation error in the small-\(t\) interval (\([\epsilon,\frac{1}{32}]\) for \(s_{\theta_1,head}\) and \([T_{trunc},\frac{1}{2}]\) for \(s_{\theta_2,tail}\)), ignoring the rest of the training set. In adjustment approaches, the model, having learned the most difficult, complicated features of the score function, is now allowed to see the full picture, minimizing (3) in its unedited, unsliced form for most optimal performance. On all experiment datasets in section 6, this approach is always able to speed up training, and sometimes results in a lower final error by escaping a local minimum that can trap a model trained with only adjustment epoches.

4 The Density Estimator↩︎

4.1 Model Design↩︎

Relying on score information, \(s_\theta(x,t)\), as input, the density estimator calculates pointwise joint densities for the selectivity estimator using the formula \[\log p_{0,\theta}(x_0)\approx\mathrm{E}_{p_{0T}(x|x_0)}\log \varphi_{\sigma^2(T)}(x)+\int_0^T n\alpha(t)dt+\frac{1}{2}\int_0^T\beta^2(t)\] \[\mathrm{E}_{p_{0t(x|x_0)}}[\left\|\nabla \log p_{0t}(x|x_0)-s_\theta(x,t)\right\|_2^2-\left\|\nabla \log p_{0t}(x|x_0)\right\|_2^2]dt,\] where we denote \[\left\|\nabla \log p_{0t}(x|x_0)-s_\theta(x,t)\right\|_2^2-\left\|\nabla \log p_{0t}(x|x_0)\right\|_2^2\] as \(g_{x_0}(x,t,s_\theta)\) for the sake of convenience. Here, \[f(x_0,\alpha,\beta,T)=\mathrm{E}_{p_{0T}(x|x_0)}\log \varphi_{\sigma^2(T)}(x)+\int_0^T n\alpha(t)dt\] can be exactly calculated with an analytic expression, so our main job would be to estimate \[\frac{1}{2}\int_0^T\beta^2(t)\mathrm{E}_{p_{0t(x|x_0)}}g_{x_0}(x,t,s_\theta)dt\] using a numerical integration scheme.
Again, due to early stopping and score blowup, our network \(s_\theta(x,t)\) is not (and cannot be) trained to produce the value of \(\nabla \log p_\delta(x)\) for all \(\delta<\epsilon\), potentially leading to wildly inaccurate representations. There are three ways to deal with this: we can evaluate the integral while skipping the \([0,\epsilon)\) interval, we can assume \(s_\theta(x,\delta)=s_\theta(x,\epsilon)\) for all \(\delta<\epsilon\), or we can define \(\tilde{s}_\theta(x,t)=s_\theta(x,t+\epsilon)\), \(\tilde{\alpha}(t)=\alpha(t+\epsilon)\), \(\tilde{\beta}=\beta(t+\epsilon)\) and calculate \[\log p_{\epsilon,\theta}(x_0)\approx f(x_0,\tilde{\alpha},\tilde{\beta},T)+\int_0^T\tilde{\beta}^2(t)\mathrm{E}_{p_{0t(x|x_0)}}g_{x_0}(x,t,\tilde{s}_\theta)dt,\] which would be an approximation of \[p_\epsilon(x)=[(\tau_{k(\epsilon)}\circ p_0)*\varphi_{\sigma^2(t)}](x)\] more solidly backed by SDE theory.
ADC chooses the third approach despite the added perturbation, both because it’s easier to analyze theoretically, and because the first two approaches seem to introduce arbitrary spikes in the approximated joint distribution function that significantly reduces the accuracy of the current selectivity estimator, designed using a predictor-corrector importance sampling algorithm. Whether the first two approaches might work for a different selectivity estimator design can be a topic of future research.
The above formula requires us to integrate \(g_{x_0}\) across both sample space (when calculating \(\mathrm{E}_{p_{0t(x|x_0)}}g_{x_0}(x,t,\tilde{s}_\theta)\)) and time. Here, time comes at a much lower dimension than space. What’s more, we generally know much more about how the norm of \(\nabla \log p_t(x)\) changes with respect to time than how it does with respect to space, as the latter greatly depends on the underlying distribution.
To account for such, we present a hybrid integration approach: for integration across time, ADC implements an adaptive stepsize Midpoint Rule Integrator, which is more suited for low-dimensional functions and can be easily modified to capitalize our knowledge of the score function’s behavior; for calculating the expectation across space, ADC implements a Quasi Monte Carlo [25] integrator much more robust against the curse of dimensionality.
Summing it up, for our chosen VP perturbation scheme where \(\alpha(t)=1\) and \(\beta(t)=\sqrt{2}\) are set to constants, ADC approximates \[\int_0^T\mathrm{E}_{p_{0t}(x|x_0)}g_{x_0}(x,t,\tilde{s}_\theta)dt\] using \[\begin{align} \sum_{i=1}^N\frac{(t_i-t_{i-1})}{2^k}\Big[\sum_{j=1}^{2^k}g_{x_0}\Big(\sigma(\tilde{t}_i)\Phi^{-1}\big(z_j+y(\tilde{t}_i)\big)+\frac{x_0}{k(\tilde{t}_i)},\tilde{t_i},\tilde{s}_\theta\Big)\Big], \end{align}\] where \(0=t_0<t_1<...<t_{N-1}<t_N=T\) is a partition of \([0,T]\), \(\tilde{t}_i=\frac{t_{i-1}+t_i}{2}\); \(\{z_j\}_{j=1}^{2^k}\) is a Sobol sequence of length \(2^k\) in the unit cube; \(\{y(\tilde{t}_i)\}_{i=1}^N\) are \(N\) random samples from \(\mathrm{U}([0,1]^d)\), used to perturb \(\{z_j\}\) so that the integrator will not focus too much on any specific region; \(\Phi\) is the distribution function of the standard normal distribution, acting pointwise on each coordinate; and the sum between \(z_j\) and \(y(\tilde{t}_i)\) is taken in \(\mathrm{T}^n\) rather than \(\mathbb{R}^n\) to prevent out-of-domain inputs for \(\Phi^{-1}\).

4.2 Integrand Choice↩︎

Based on the data prediction model, we further propose modifying the expression of \(g_{x_0}(x,t,\tilde{s}_\theta)\) for \(t>T_{trunc}\) in a way that drastically reduces integration error and allows for the use of much larger time-steps.
Our idea is based on the following observations. First, for large \(t\), \(s_\theta(x,t)\) would quickly approach \(-\frac{x}{\sigma^2(t)}\). Second, the data prediction model is already using \(\frac{v_\theta(x,t)}{\sigma^2(t)k(t)}\) to approximate the residual term \(\nabla \log p_t(x)+\frac{x}{\sigma^2(t)}\) rather than \(\nabla \log p_t(x)\) itself. Most importantly, \[\mathrm{E}_{p_{0t}(x|x_0)}\big[\left\|\nabla \log p_{0t}(x|x_0)+\big(x/\sigma^2(t)\big)\right\|_2^2-\left\|\nabla \log p_{0t}(x|x_0)\right\|_2^2\big]\] is exactly calculated as \[\frac{\beta^2(t)(n\sigma^2(t)k^2(t)+\left\|x_0\right\|_2^2)}{\sigma^4(t)k^2(t)},\] which can be integrated analytically on any interval whenever \(\alpha(t)\) and \(\beta(t)\) are constants.
With the above in mind, for all \(t>T_{trunc}\), ADC decomposes \(g_{x_0}(x,t,\tilde{s}_\theta)\) into the terms \[\frac{1}{\sigma^4(t)k^2(t)}\langle\tilde{v}_{\theta}(x,t),\tilde{v}_{\theta}(x,t)-2x_0\rangle-\langle\frac{x}{\sigma^2(t)},\frac{x}{\sigma^2(t)}-\frac{2x_0}{\sigma^2(t)k(t)}\rangle,\] leaves the second term for exact calculation, and only approximates the first term, denoted from now on as \(g_{res,x_0}(x,t,\tilde{v}_\theta)\), using hybrid numerical integration.
This modification is meaningful in that as \(t\) gets larger, \(g_{x_0}(x,t,\tilde{s}_\theta)\), the original integrand, would approach a nonzero value that varies notably with \(x\), but \(g_{res,x_0}(x,t,\tilde{v}_\theta)\) converges to zero at a rate no slower than \(\mathcal{O}(\frac{1}{\sigma^4(t)k^2(t)})\), i.e., \(\mathcal{O}(e^{-2t})\) under the VP scheme. As a result, ADC is allowed to implement much larger time steps for even moderately large \(t\), significantly boosting its efficiency.

4.3 Parameter Choice↩︎

We discuss our rationale for choosing two parameters: \(T\), the integration upper bound, currently set to 3; and \(\{t_i\}_{i=0}^N\), the timestep scheme, currently chosen via the pseudocode in Algorithm 1.
To choose \(T\), we consider the following theorem that provides a lower bound for the speed at which \(p_t\) converges to a known Gaussian distribution \(\mathcal{N}(0,\sigma^2(t)I_d)\).

The above theorem tells us that under the condition \(\mathrm{E}_{p_0(x)}x=0\), predicting \(\nabla \log p_t(x)\) using the score function of a normal distribution results in an error term that decays at a speed of roughly \(\mathcal{O}(\frac{1}{\sigma^4(t)k^2(t)})\) for small \(t\), and roughly \(\mathcal{O}(\frac{1}{\sigma^6(t)k^4(t)})\) for large \(t\), with the transition guaranteed to happen when \[\frac{2\left\|x_i\right\|_2^2+\mathrm{E}_{p_0(x)}\left\|x\right\|_2^2}{2\sigma^2(t)k(t)}\leq \log 2\] holds for all \(\{x_i\}_{i=1}^M\) in the database. Meanwhile, one can easily deduce from (5) that \[D_{KL}(p_T(x),\mathcal{N}(0,\sigma^2(T)))=\int_T^\infty\mathrm{E}_{p_t(x)}\left\|\nabla\log p_t(x)+\frac{x}{\sigma^2(t)}\right\|_2^2dt,\] so setting \(T\) slightly beyond this transition point might be the most cost-effective approach.
In the case of ADC, our normalization process ensures \(\mathrm{E}_{p_0(x)}=0\), and each attribute’s \(\max\) and \(\min\) value differ by 3.2. As a result, under the VP scheme, setting \(T=3\) would suffice for all datasets of moderate dimensionality.
To choose the timestep scheme, we first perform an error analysis on our hybrid integrator. For an arbitrary timestep partition scheme \(\{t_i\}_{i=0}^N\), define the time-related discretization error \(\delta_{i,time}(x_0)\) as \[\int_{t_i}^{t_{i+1}}[\mathrm{E}_{p_{0t}(x|x_0)}g_{x_0}]dt-\Delta t_i\mathrm{E}_{p_{0\tilde{t}}(x|x_0)}g_{x_0}(x,\tilde{t}_i,\tilde{s}_\theta),\] and the space-related discretization error \(\delta_{t,space}(x_0,y)\) as \[\mathrm{E}_{p_{0t}(x|x_0)}g_{x_0}(x,t,\tilde{s}_\theta)-\frac{1}{2^k}\sum_{j=1}^{2^k}g_{x_0}\big(\sigma(t)\Phi^{-1}\big(z_j+y\big)+\frac{x_0}{k(t)},t,\tilde{s}_\theta\big).\] Given the curse of dimensionality and that the volatile \(g(x,t,s_\theta)\) was already "smoothed up" somewhat upon calculating \(\mathrm{E}_{p_{0t}(x|x_0)}g_{x_0}\), one can naturally expect \(\{\delta_{i,time}(x_0)\}_{i=1}^N\) to be significantly smaller than \(\{\delta_{\tilde{t}_i,space}(x_0,y(\tilde{t_i}))\}_{i=1}^N\) with a dense enough timestep scheme, in which case we have the following theorem.

One issue yet unresolved is that \(\mathrm{E}_{x_0\sim p_0}[\mathrm{Var}_{y}\delta_{\tilde{t}_i,space}(x_0,y)]\) proves very hard to estimate. As a result, the current version of ADC (somewhat boldly) deems the term proportionate to \[\mathrm{E_{x_0\sim p_0}\mathrm{E}_{p_{0t(x|x_0)}}}\big(g_{x_0}(x,t,\tilde{s}_\theta)\big)^2\] when \(t<T_{trunc}\) and \[\mathrm{E_{x_0\sim p_0}\mathrm{E}_{p_{0t(x|x_0)}}}\big(g_{res,x_0}(x,t,\tilde{s}_\theta)\big)^2\] when \(t>T_{trunc}\), leaving the task of more accurately estimating \(\mathrm{E}_{x_0\sim p_0}[\mathrm{Var}_{y}\delta_{\tilde{t}_i,space}(x_0,y)]\) for upcoming research.
The decay rate of the second term is estimated to be no slower than \(\mathcal{O}(\frac{1}{\sigma^8(t)k^4(t)})\) when \(t\to\infty\), as we already discussed. For the first term, Theorem 1 proves that for all \(x_0\) such that \(p_0(x_0)\neq 0\), \(\mathrm{E}_{p_{0t}(x|x_0)}\left\|\nabla\log p_{t+\epsilon}(x)\right\|_2^4\) can be upper bounded by \[8\mathrm{E}_{p_{0t}(x|x_0)}\left\|\frac{x-x_0}{\sigma^2(t+\epsilon)}\right\|_2^4+8\mathrm{E}_{p_{0t}(x|x_0)}\left\|\delta_{\sigma^2}(x,t+\epsilon)\right\|_2^4,\] which, using the bounds on \(\delta_{\sigma^2}(x,t)\), can be roughly estimated as proportionate to \(\frac{1}{\sigma^4(t+\epsilon)}\). Therefore, Cauchy’s Inequality proves the growth rate of \[\mathrm{E}_{p_{0t}(x|x_0)}(\langle\nabla\log p_{t+\epsilon}(x),\nabla\log p_{t+\epsilon}(x)-2\nabla\log p_{0t}(x|x_0)\rangle)^2\] to be no faster than \(\mathcal{O}(\frac{1}{\sigma^2(t)\sigma^2(t+\epsilon)})\) as \(t\to 0\), thus justifying our pseudocode implementation.

5 The Selectivity Estimator↩︎

5.1 Predictor-Corrector Integration: Why GMM↩︎

While the score estimator and density estimator are effective in calculating the probability density function’s value at any given point, we still need to integrate this function over the entire queried region in order to estimate the selectivity of ranged queries. ADC does this by utilizing the importance sampling Monte Carlo formula \[\int_V p(x)dx=\int_Vq(x)dx\cdot\mathrm{E}_{x\sim q(x|x\in V)}\frac{p(x)}{q(x)}\] that Naru [4] and DQM-D[12]’s selectivity estimator are both based on, which, when the expectation on the right is approximated using i.i.d. sampling, becomes an unbiased estimator with variance proportionate to \[\Big(\int_Vq(x)dx\Big)^2\mathrm{Var}_{x\sim g(x| V)}\Big(\frac{p(x)}{q(x)}\Big).\] To enhance speed and precision, the selectivity estimator of ADC must construct a predictor that can be rapidly integrated, rapidly sampled from, and resembles the corrector \(p_{\epsilon,\theta}\) as much as possible. However, the predictor constructed by Naru (via progressive sampling) and DQM-D (via VEGAS [13]) require multiple phases of sampling. This is undesirable for us, considering that evaluating the density function at one point already requires a not-so-easy numerical integration process.
Gaussian Mixture Models (GMMs) attempt to reconstruct a distribution \(q\) using the expression \[\begin{align} q=\sum_{i=1}^Nw_i\mathcal{N}(y_i,H_i), \end{align}\] and ADC’s choice to train one as its predictor is born out of the very same considerations listed above.
In terms of speed, the advantage of GMM manifests in that for \(\{H_i\}_{i=1}^N\) diagonal and \(V=\{(x_1,..,x_d)|a_j\leq x_j\leq b_j, \forall j\}\) a hyper-rectangle in \(\mathbb{R}^d\), \(\int_Vq(x)dx\) can be calculated as \[\int_Vq(x)dx=\sum_{i=1}^Nw_i\prod_{j=1}^d\left[\Phi(h_{jj}(b_j-y_j))-\Phi(h_{jj}(y_j-a_j))\right],\] and conditional samples \(\{x_j\sim q(x|V)\}\) can be drawn using \(d+1\) random numbers \(z_i\sim U[0,1]\) per sample, with \(z_{n+1}\) determining which kernel \(\mathcal{N}(y_j,H_j)\) \(x\) belongs to, and \(z_i\) determining the \(i\)th coordinate of \(x\) to be \[x_i=\Phi^{-1}\big(z_i\Phi(h_{ii}(b_i-y_{j,i}))+(1-z_i)\Phi(h_{ii}(y_{j,i}-a_i))\big).\] As a result, when \(V\) is a hyperrectangle, sampling and integration can both be performed with a time cost independent of \(V\)’s diameter and comparable to that of calculating a pointwise density value.
In terms of similarity to \(p\), another advantage of GMM is that the underlying distribution \(p_\epsilon\), which \(p_{\epsilon,\theta}\) seeks to approximate, always have a KDE-type expression \[p_\epsilon (x)=\frac{1}{M}\sum_{i=1}^M\mathcal{N}(k(\epsilon)^{-1}(x_i),\sigma^2(\epsilon)),\] as in the database case, \(p_0\) is simply the average of delta functions placed at sample points. Given their unique shape and smoothness features, "Gaussian mixtures work best at matching local features of Gaussian mixtures" is a very reasonable assumption.
Instead of a query-based training approach, ADC trains its GMM predictor to maximize the log-likelihood function \[\begin{align} \frac{1}{M}\sum_{j=1}^M\log q(x_j). \end{align}\] using the EM algorithm [26], an elegant GMM-exclusive optimizer that trains GMMs by iteratively alternating between estimating the posterior responsibilities of mixture components (E-step) and updating the model parameters to maximize the expected complete-data log-likelihood (M-step). Compared to one-size-fit-all optimizers like Stochastic Gradient Descent (SGD), EM is very robust and efficient, as it requires no hyperparameters, guarantees monotone increase of likelihood, and always converges to a local maximum value.
ADC further attempts to escape from local maximums, something the classical EM algorithm is incapable of, by resampling kernels with extremely low weights. That is, after every 10 EM iterations, ADC deletes all Gaussian kernels with weights below \(1/500N\), and replace them with new kernels of weight \(1/500N\), covariance \(k^2I_d\), and mean values resampled from \(\{x_j\}_{j=1}^M\) via a weighted sampler that sets the weight of a data point \(x_i\) inversely proportionate to the current value of \(q(x_i)\). If training time is sufficient, it is also sometimes preferable to train multiple GMMs at once, and then select the model that achieves the highest likelihood.

5.2 Neutralizing Early Stopping Error↩︎

We briefly discuss two ways of using our GMM predictor-corrector scheme to cancel out the error caused by early stopping, shown by experiments to be crucial for ADC’s optimal performance.
Our first method is based on the following observations: First, \(q\), when trained to maximize the function (8), is in fact trying to approximate not \(p_\epsilon\), but \(p_0=\frac{1}{M}\sum_{i=1}^M\delta(x_i)\) itself. Second, perturbing GMMs using a diffusion process is extremely easy: setting \(q_0(x)\) to be the \(q(x)\) defined in (7), we would simply have \[q_t=\sum_{i=1}^Nw_i\mathcal{N}\Big(\frac{y_i}{k(t)},\frac{H_i}{k(t)}+\sigma^2(t)I_d\Big).\] With this in mind, ADC’s selectivity estimator assumes that \[\frac{\int_Vp_{0,\theta}(x)dx}{\int_Vp_{\epsilon,\theta}(x)dx}\approx\frac{\int_Vq_0(x)dx}{\int_Vq_\epsilon(x)dx},\] and then corrects for the early stopping error by actually calculating \[\int_Vp_{0,\theta}(x)dx\approx\int_Vq_0(x)dx\cdot\mathrm{E}_{x\sim q_\epsilon(x|V)}\Big(\frac{p_{\epsilon,\theta}(x)}{q_\epsilon(x)}\Big).\] Our second method is based on the observation that both GMMs and diffusion models work by "bleeding" portions of density from the typical set onto neighboring regions. When one or more pairs of attributes \((A_i,A_j)\) exhibit near-functional dependency, hence the typical set’s projection onto \((A_i,A_j)\) becomes a 1-dimensional curve, such a practice significantly distorts the local probability distribution, leading to large errors for extreme-case queries.
ADC mitigates this problem by using Bayesnet [19] to handle functionally dependent attribute pairs. In case where knowing the value of an attribute \(A_i\) can help narrow down the range of another attribute \(A_j\) to \(k\) (chosen to be 0.05 in the current model) times the original, ADC, instead of training a GMM and a diffusion model on all attributes, would instead work by training:

• a diffusion model \(p_\epsilon\) and a GMM \(q\) modeling the distribution of all attributes excluding \(A_j\)

• a Bayesnet histogram modeling the conditional probability density \(p(A_j|A_i)\) as a function piecewise constant with respect to \(A_i\) and \(A_j\).

Now, given a query asking for all points inside \(V=\prod_{k=1}^d[a_k,b_k]\), ADC predicts its selectivity using \[\sum_{m=1}^Nw_m\left(\prod_{k\neq{i,j}}\int_{a_k}^{b_k}\varphi_{H_{m,kk}}(x-y_{m,k})dx\cdot g(m,a_j,b_j)\right),\] with \(g(m,a_j,b_j)\) defined as \[\int_{a_i}^{b_i}\varphi_{H_{m,ii}}(x-y_{m,i})p(y_j\in[a_j,b_j]|x_i)dx_i,\] then corrects its prediction using the weighed average of \(\{\frac{p_\epsilon(z_k)}{q_\epsilon(z_k)}\}_{k=1}^n\), with \(\{z_k\}_{k=1}^n\) i.i.d. samples from \(q_\epsilon\) and the weight of each sample point set proportionate to \[\{p(x_j\in [a_j,b_j]|z_{k,i})\}_{k=1}^n.\] Because \(p\), constructed using histograms, is constant with respect to \(x_i\) in each of the intervals \(\{[x_{i,p},x_{i,p+1}]\}_{p=0}^{P-1}\), calculation of the above expressions can be greatly sped up by calculating and caching, for the dimension \(i\), each slice \([x_{i,p},x_{i,p+1}]\), and each Gaussian kernel \(m\), \[\int_{x_{i,p}}^{x_{i,p+1}}\varphi_{H_{m,ii}}(z_{m,i}-x_i)dx_i\] once upon initializing the model; as well as calculating \[\{p(y_j\in[p_j,q_j]|x_i\in[x_{i,p},x_{i,p+1}])\}_{p=0}^{P-1}\] once upon receiving a query and using the values for every Gaussian kernel.
For our benchmark datasets, the above Bayesnet implementation worked very well on power, successfully cutting down the max Q-error of queries from 100 to 28, as the range of the 0th attribute can be narrowed down to an average of 4.5% its original by knowing the 4th attribute’s value.

5.3 Introducing ADC+↩︎

We briefly discuss ADC+, another GMM-based addition aimed at enhancing ADC’s speed and accuracy on large-volume, high-selectivity queries.
In numerical experiments, we observe that the GMM predictor itself can handle such "easy" queries quite competently, often resulting in a level of precision that rivals or even exceeds the entire prediction-sampling-correction algorithm. Thus, we further upgrade ADC by training a decision tree classifier \(\mathrm{T}: \mathbb{R}^2\to \{0,1\}\) and implementing the pseudocode in Algorithm 2 (for conciseness, we depict the general case where Bayesnet is not used; in case it is, the algorithm can be easily modified by replacing the above formulas with those in 5.2).
ADC+ trains \(T\) to minimize the square of log Q-error. That is, given a set of sample queries, we use gini impurity as our loss function and set the weight of each query to \(|(\log Q_{ADC})^2-(\log Q_{GMM})^2|\), where \(Q_{ADC}\) and \(Q_{GMM}\) refers to the Q-error on that query when estimated using the entire ADC algorithm or only the GMM predictor, respectively. To prevent overfitting and reduce storage costs, the decision tree is designed to always have a maximum depth of 4.
Considering that the advantages of ADC+ are more easily verified experimentally instead of theoretically, more information regarding our motives behind creating ADC+, as well as the performance gains resulting from it, shall be provided in 6.3.

6 Experiments↩︎

6.1 Experimental Setup↩︎

6.1.1 Datasets↩︎

We test the performance of ADC on two real world and one synthetic dataset, including:
Forest[5]: Forest coverage type dataset with 581012 rows, used as a benchmark by multiple influential works on cardinality estimation [2] [17] [6] [8]. Following the practice of [2], we keep the first 10 continuous attributes for evaluation.
Power[5]: Dataset with 2.1 million rows, recording power consumption information of a French household at a 1 minute sampling rate, used alongside forest in all four essays above. Following the practice of [2], we keep the 7 continuous measurement attributes excluding date and time, and convert "?" to -1 when dealing with the 1.25% of rows consisting entirely of missing values.
Modulo: A synthetic dataset designed according to the philosophy that lead to our creation of ADC: treating all attributes as a single, unified entity is at times necessary for producing decent estimates. The dataset consists of 4 million rows and five attributes \(A\),\(B\),\(C\),\(D\), and \(E\), where \(A\), \(B\), and \(D\) are independent integers taking the values 0 to 1999 with equal probability; \[C=(A+B+\epsilon_1)\bmod 2000, E=(A+D+\epsilon_2)\bmod 2000,\] and \(\epsilon_1\), \(\epsilon_2\) are independent integers taking the values 0 to 199 with equal probability. A property of modulo is that all five attributes are pairwise independent from one another, making it potentially challenging for estimators that work by discerning pairwise correlations.

6.1.2 Workloads↩︎

For each dataset, we create 100,000 queries for training query-driven methods, 10,000 queries for validation, and 10,000 queries for testing using the program provided by [2]. That is, for a \(d\)-dimensional dataset, we design queries asking for all points lying within a hyperrectagle, where:

• The hyperrectangle works as a filter for \(k\) attributes and spans the entire range of all others, with \(k\) a random integer uniformly selected between 1 and \(d\).

• The center of the hyperrectangle has a 90% chance of being randomly sampled from the dataset (data-entric), and a 10% chance of being randomly sampled from a uniform distribution spanning the range of each attribute (random-centric)

• For each attribute with a non-trivial range bound, the query range width on that attribute has a 50% chance of being sampled from \(\mathrm{U}([0,R_i])\) and a 50% chance of being sampled from \(\mathrm{Exp(\frac{10}{R_i})}\), with \(R_i\) denoting the difference between \(\max\) and \(\min\) values of that attribute.

6.1.3 Comparative Techniques↩︎

We briefly introduce versions of our model and the models we use for comparison, which include all learned estimators tested in [2]. For the sake of giving each model a similar size budget while preventing poor results due to bad parameter choices, comparative techniques are trained using the same hyperparameter settings the author of [2] chose for the dataset forest, putting the parameter size of each model at around 0.5MB per dataset.
ADC: Our estimator ADC in its raw form, always going through the prediction-sampling-correction procedure for every query.
ADC+: ADC in its most optimal form, after implementing the changes proposed in 5.3.
Pred: The performance of ADC’s GMM predictor as a standalone estimator, which we report to to quantify accuracy improvements gained by the diffusion model corrector, i.e. the score estimator and density estimator.
MSCN [6]: a query-driven regression model which processes the feature vector of queries using a multi-set convolutional network.
Lw-Tree and LW-NN [8]: two lightweight query-driven regression models that process the feature vector of queries using either a gradient boosted tree (LW-Tree) or a neural network (LW-NN), and enhance their precision with the help of CE (correlation based) features extracted from simple heuristic estimators.
DeepDB [7]: a data-driven joint distribution model that reconstructs data distributions using sum-product networks.
Naru [4]: a data-driven joint distribution model that reconstructs data distributions by training an autoregressive model and answers ranged queries using sequential sampling, shown by [2] to be the most robust and accurate learned model on most test datasets.

6.1.4 Experiment Device and Accuracy Metric↩︎

We run all our experiments on a dual-socket Intel Xeon Gold 6140 system with 36 physical cores (72 hardware threads) and a base clock frequency of 2.30 GHz. We choose Q-error, defined as \[\max\{\frac{Card_{real}}{Card_{est}},\frac{Card_{est}}{Card_{real}}\}\] to be our accuracy metric, with zero values for \(Card_{real}\) and \(Card_{est}\) both set to 1 to prevent division by zero. We report \(\max\), 99th percentile, 95th percentile, median, and geometric mean Q-error as our evaluation benchmarks in accordance with [2].
We note that the performance of LW-Tree/LW-NN on forest and MSCN, Naru, DeepDB on power are not as decent in our experiments as they were in [2]’s, possibly due to different Python library versions our smaller model size budget for the power dataset⁴. While we record in Table 1 the results from our experiments for methodological consistency, we also encourage readers to look up the original figures in [2] for a more comprehensive comparison.

6.2 Comparative Experiment Evaluation↩︎

6.2.1 Accuracy Evaluation↩︎

The results for all estimators are listed in Table 1(a) for forest, power and Table 1(b) for modulo, which we shall now discuss.
We begin by discussing the real-world datasets. On the most important benchmarks GM, 95th, and 99th Q-error, ADC+ and Naru consistently demonstrate top-tier performance, rivaling each other while significantly outperforming every other comparative technique. Meanwhile, ADC performs almost identically to ADC+ on 95th and 99th Q-error, but falls somewhat behind on GM Q-error due to its higher median Q-error values. As Naru, the only estimator which sometimes outperform ADC+, comes with the downside of much higher latency, such experimental results are enough to earn ADC+ a place among the best cardinality estimators.
Regarding max Q-error, the performance of all six estimators seem to be highly volatile and significantly influenced by random factors: for example, on forest, DeepDB performed notably better in our experiment than it did in [2], but its max Q-error increased from 1293 to 3856; on power, LW-NN performed roughly the same as in [2], but its max Q-error dropped from 40000 to 22000; meanwhile, while ADC’s largest Q-error was 551 on forest and 267 on power, it’s second-largest Q-error dropped sharply to 137 and 46.9, respectively, meaning that its max Q-error might change drastically on a slightly different workload. That said, ADC and ADC+ are still capable of outperforming every model except Naru and MSCN on max Q-error, proving ADC to be a relatively robust estimator.
For the dataset modulo, ADC’s advantage scenario, the complex, multi-attribute correlations caused Naru, DEEPDB, and MSCN to significantly drop in accuracy, while LW-Tree’s performance exhibited less drastic declines. LW-NN is the only estimator to exhibit performance comparable to ADC on this dataset, almost catching up with ADC+ on 95th and 99th Q-error but still falling significantly behind ADC+ on the benchmarks GM and max. Meanwhile, the GMM predictor of ADC, a close relative of the KDE algorithm, only exhibits barely satisfactory performance on 95th and 99th Q-error, while Bayesnet implemented via the Chow-Liu algorithm, which only captures pairwise correlations, can be proven incapable of learning anything more than single-attribute histograms, showing that this performance decline is not unique to learned cardinality estimators.
Except for max Q-error on power, ADC consistently outperforms its standalone GMM predictor on GM, 95th, 99th, and max Q-error, but the magnitude of this improvement varies greatly across datasets: ADC reduces the distance of 99th Q-error to 1 by roughly one half on modulo, one fourth on forest, and one sixth on power. Thus, except for scenarios with intense latency requirements, using diffusion models for cardinality estimation will always pay off.

6.2.2 Latency and Storage Cost↩︎

To evaluate the practicality of ADC+ in actual user scenarios, we discuss the inference latency and storage requirements of ADC+ in comparison with other learned models.
In terms of storage space, ADC+ is at an advantage compared to all other models, with parameters taking up 0.34MB of space for forest and less than 0.3MB for power and modulo. To be more precise, the network \(s_{\theta_1,head}\) usually takes up around 150 KB of storage, \(s_{\theta_2,tail}\) around 50 KB, the GMM predictor usually takes up around \(10d\) KB with \(d\) the dataset’s dimensionality, and other artifacts, including the decision tree, some database information, and a histogram we use for single-attribute queries usually take up around 20 KB. Meanwhile, all other models in our experiments have at least 0.42MB of trainable parameters (except Naru on power), and the actual number often goes further beyond that. We note that other learned models can be made more lightweight if needed, though often at a cost. For example, in their original essay, LW-Tree and LW-NN [8] conducted their experiments using models as small as 16KB but, as a result, measured significantly higher Q-errors on forest and power than we did.
Regarding per-query latency, the most important feature besides accuracy, ADC+ is much slower than MSCN, LW-Tree and LW-NN, which answer a query within 0.5 ms (within 0.2ms for LW-Tree), and somewhat slower than DeepDB, but is also \(2\times\) as fast as Naru on forest and power while more than \(4\times\) as fast as Naru on modulo. We find that ADC+’s slower speed on forest might be due to two factors: the heavier network we implemented due to the higher dimensionality and looser correlations, and the fact that ADC’s GMM predictor becomes less accurate as dimensionality rises, producing a satisfactory estimation by itself (hence allowing ADC+ to skip the sampling and correction phases) only on very rare occasions. Given that ADC+ is almost entirely written in Pytorch, it might also be possible to further speed up ADC+ with the help of GPUs which, sadly, we do not have access to.

6.3 Reducing Median Q-error: Why ADC+↩︎

Our greatest motivation for creating ADC+ is that ADC, in its raw form, performs somewhat awkwardly on median Q-error, nearly always losing to its GMM predictor. In fact, on power and modulo, ADC ironically increased the distance between median Q-error and 1 by five times, ranking embarrassingly near bottom place on this benchmark despite its decent 95th, 99th, and max Q-error performance.
Why is this? We suspect that it’s caused by the fact that unlike other estimators, there are no "very easy" queries for ADC, and even a query that only asks about two independent attributes or one whose queried region spans 90% of sample space must undergo the same prediction-sampling-correction procedure. In this case, ADC’s predictor-corrector integration scheme often finds itself adding non-negligible variance to correct an already tiny bias, a choice that is increasingly likely to backfire as the GMM predictor grows more competent.
To test our hypothesis, we create a scatterplot documenting the query volume and GMM-predicted selectivity of 100,000 training queries for each of the two real-world datasets, with blue points depicting queries where the GMM predictor performs better than ADC and red points depicting the opposite. Indeed, we can see from the experimental results documented in Figure 3 that while red points appear more frequently than blue in most regions of the graph, each graph’s upper-right corner always contains a very dense cluster that is almost exclusively blue.
What’s even better, separating a large proportion of blue points from red, at least those on the upper-right corner, can be done by drawing simple, clear-cut boundary lines: simple enough to be handled by a decision tree, which is exactly what ADC+ does. Compared to neural networks, a decision tree of depth 4 only takes up less than 3KB and classifies a query in less than 20 microseconds, allowing it to be implemented almost for free.
One might worry whether wrongly classifying points that actually benefit from predictor-corrector integration would lead to worse 95th and 99th percentile Q-error results. Fortunately, experiments show that ADC+’s effects on these benchmarks, though small, is almost never negative. On median Q-error, meanwhile, ADC+ acts as the game changer, successfully beating MSCN, LW-Tree, and LW-NN due to roughly halving the distance between median Q-error and 1 on all three datasets. What’s even better, ADC+ is also able to cut down the per-query latency of ADC by roughly 25%, because the benefits gained from potentially avoiding two lengthy steps of numerical integration often far outweighs the negligible classification costs.

7 Conclusion↩︎

The ability of diffusion models to reconstruct complex, high dimensional distributions makes them accurate cardinality estimators. We have proven this by constructing ADC, whose accuracy exceeds many of its famous peers on both real-world and synthetic datasets. Downsized using newfound theorems regarding the score function’s behavior, ADC also takes up less than 400KB of storage, and keeps its latency within 10ms on two out of three datasets tested.
ADC’s current capabilities are still somewhat limited by its inability to handle categorical attributes, an issue we seek to resolve in future studies. We are also curious whether some theorems, network structures, and modeling approaches we proposed upon designing ADC, such as switching from noise to data prediction as \(t\) gets larger, might be beneficial for training diffusion models in general, despite influences of the manifold hypothesis and the much higher (3000+) dimensionality.
Overall, ADC provides a good starting point for exploring the use of diffusion models in cardinality estimation, a rather unconventional crossover that can produce very promising fruits.

8 Mathematical Notations Used in the Essay↩︎

9 Proof of Theorems and Claims↩︎

::: theorem 3.1 Theorem 3.1 1. Let \(p_0\) be the weighted average of several positive distributions \(\{p_i\}_{i=1}^n\), such that \[\{V_i\}_{i=1}^n:=\{x|p_i(x)>0\}_{i=1}^n\] are open domains in \(\mathbb{R}^d\) intersecting only at their boundaries, and all \(p_i\) are Lipschitz with constant \(k\) inside \(\bar{V_i}\). Denote \(V_0=\mathbb{R}^d\textbackslash\cup_{i=1}^n \bar{V_i}\) for convenience, and let \(\{p_t|t\in(0,T]\}\) be the family of distributions derived by perturbing \(p_0\) according to the VE scheme (\(\sigma^2(t)=t, k(t)=1\)), then:

1. \(\forall x_0\in V_i\) such that \(B(x,r)\subseteq V_i\), if \(\nabla\log p_i(x)\) exists and is Lipschitz continuous in a neighborhood of \(x_0\) with constant \(h\), then \[\nabla \log p_t(x_0)=\nabla \log p_0(x_0)+\delta_1(x_0,t),\] with \(\delta_1(x_0,t)\) an infinitesimal term of order \(O(\sigma(t))\).

2. \(\forall x_0 \in \partial V_i\cap\partial V_j\) with at most one of \(i\) and \(j\) being zero, if \(B(x_0,r)\subseteq V_i\cup V_j\), \(p_{0,i}(x_0)+p_{0,j}(x_0)>0\), and inside \(B(x_0,r)\), \[\mathrm{d}\left(y,T_{x_0}(\partial V_i\cap \partial V_j)\right)\leq h\left\|y-x_0\right\|_2^2\] holds for some constant \(h\) and all \(y\in\partial V_i\cap\partial V_j\), then \(\forall \lambda\), \(\nabla \log p_t(x_0+\lambda\sqrt{t}\vec{n})\) would equal \[\frac{e^{-\frac{\lambda^2}{2}}(p_{0,i}(x_0)-p_{0,j}(x_0))\vec{n}}{2\sigma(t)\left(\Phi(\lambda)p_{0,i}(x_0)+(1-\Phi(\lambda)p_{0,j}(x_0))\right)}+\delta_\sigma(x_0,\lambda,t),\] with \(\vec{n}\) the unit vector normal to \(T_{x_0}(\partial V_i\cap\partial V_j)\), pointing from \(V_j\) to \(V_i\), and \(\delta_\sigma(x_0,\lambda,t)\) a residual term of order \(O(1)\).

3. \(\forall x_0\in V_0\), suppose \(y_0\in \bar{V_i}\) is the point closest to \(x_0\) such that \(y_0\in\cup_{i=1}^k\bar{V_i}\). If \(P(x\in B_r(y)|x\sim p_0)>M_1r^{k_1}, \forall 0<r<R\) for some constants \(M_1, k_1\), while \[\left\|z-x_0\right\|_2^2-\left\|y_0-x_0\right\|_2^2\geq \min\{M_2,k_2\left\|z-y_0\right\|_2^2\}\] for all \(z\in\cup_{i=1}^kV_i\) and some constants \(M_2,k_2\), then \[\nabla\log p_t(x_0)=\frac{y_0-x_0}{\sigma^2(t)}+\delta_{\sigma^2}(x_0,t),\] with \(\delta_{\sigma^2}(x_0,t)\) a residual term of order \(o(\frac{1}{\sigma(t)^{1+\epsilon}})\) for all \(\epsilon>0\). :::

::: theorem 4.1 Theorem 4.1 1. Let \(\tilde{x}\) be a random variable taking values \(\{x_i\}_{1\leq i\leq n}\) with equal probability (as is always the case when we construct \(p_t\) from i.i.d. sample points), such that \(\frac{1}{n}\sum_{i=1}^nx_i=0\). Take \(p_t(x)\) to be the marginal distribution of \(x(t)\) when \(x(0)=\tilde{x}\) is perturbed by the diffusion SDE \[dx=-\alpha(t)xdt+\beta(t)dw\] from time 0 to time \(t\). Then, we would have

1. \(\mathrm{E}_{p_t(x)}\left\|\nabla\log p_t(x)+\frac{x}{\sigma^2(t)}\right\|_2^2\leq \frac{f(t,\tilde{x})}{\sigma^4(t)k^2(t)}.\)

2. \(\mathrm{E}_{p_t(x)}\left\|\nabla \log p_t(x)+\frac{x}{\sigma^2(t)}\right\|_2^2\leq \frac{f(t,\tilde{x})g(t,\tilde{x})}{\sigma^4(t)k^2(t)}.\)

where \[f(t,\tilde{x})=\frac{1}{n}\sum_{i=1}^n\left\|x_i\right\|_2^2\] and \[g(t,\tilde{x})=\frac{1}{n}\sum_{i=1}^n(e^{\frac{2\left\|x_i\right\|_2^2+f(t,\tilde{x})}{2\sigma^2(t)k^2(t)}}-1).\] :::

::: theorem 4.2 Theorem 4.2 1. Define \(\mathrm{Var}(x_0)\) to be the error term’s variance upon approximating \(\log p_{\epsilon,\theta}(x_0)\) with formula (6). Assuming that

1. The partition is dense enough that for the same \(y\), fluctuations of \(\Delta t_i\delta_{t,space}(x_0,y)\) for different \(t\) in \([t_{i-1},t_{i}]\) is negligible compared to the total integration error

2. The error terms \(\{\delta_{i,time}(x_0)\}_{i=1}^N\) are negligible compared to \(\{\delta_{\tilde{t}_i,space}(x_0,y(\tilde{t_i}))\}_{i=1}^N\)

3. The sequence \(\{y(\tilde{t}_i)\}_{i=1}^N\) consists of i.i.d. samples from \(\mathrm{U}([0,1]^d)\)

hold for all \(x_0\), \(\mathrm{E}_{x_0\sim p_0}\mathrm{Var}(x_0)\) would be minimized under a timestep scheme where \(t_i-t_{i-1}\) is chosen inversely proportionate to \[\big(\mathrm{E}_{x_0\sim p_0}[\mathrm{Var}_{y}\delta_{\tilde{t}_i,space}(x_0,y)]\big)^{\frac{1}{2}}.\] :::

References↩︎