Fast Algorithms for the Multi-dimensional Jacobi Polynomial Transform

We use the well-known observation that the solutions of Jacobi's differential equation can be represented via non-oscillatory phase and amplitude functions to develop a fast algorithm for computing multi-dimensional Jacobi polynomial transforms. More explicitly, it follows from this observation that the matrix corresponding to the discrete Jacobi transform is the Hadamard product of a numerically low-rank matrix and a multi-dimensional discrete Fourier transform (DFT) matrix. The application of the Hadamard product can be carried out via $O(1)$ fast Fourier transforms (FFTs), resulting in a nearly optimal algorithm to compute the multidimensional Jacobi polynomial transform.

Introduction
The one-dimensional forward discrete Jacobi transform consists of evaluating an expansion of the form where P (a,b) ν denotes the order-ν Jacobi polynomial of the first kind corresponding to parameters a and b, at a collection of points {x i } i=1,...,n ⊂ (−1, 1). The one-dimensional inverse discrete Jacobi transform is the process of computing the coefficients {α ν } i=1,...,n in an expansion of the form (1) given its values at a collection of distinct points {x i } i=1,...,n ⊂ (−1, 1). For the sake of brevity, we will generally drop the adjective "discrete" and simply use the terms "forward Jacobi transform" and "inverse Jacobi transform" when referring to these operations. Often, the points {x i } i=1,...,n ⊂ (−1, 1) are the nodes of the n-point Gauss-Jacobi quadrature rule p(x ν )ω ν (2) that is exact when p is a polynomial of degree less than or equal to 2n − 1. However, it is useful to consider more general sets of points as well. In the case in which the points {x i } i=1,...,n ⊂ are the nodes of the Gauss-Jacobi quadrature rule, we say the corresponding transform and inverse transform are uniform; otherwise, we describe them as nonuniform. To form an orthonormal matrix to represent the transform for numerical purpose, we will rescale the transform in (1) with quadrature weights later. Many methods for rapidly applying the Jacobi transform and various special cases of the Jacobi transform have been developed. Examples include the algorithms of [1,2,3,4] for applying the Legendre transform, and those of [5] for forming and evaluating expansions for Jacobi polynomials in general. See also, [6] and its references for extensive information on numerical algorithms for forming and manipulating Chebyshev expansions. Almost all such algorithms can be placed into one of two categories. Algorithms in the first category, such as [7,8,9,10,11], make use of the structure of the connection matrices which take the coefficients in the expansion of a function in terms of one set of Jacobi polynomials to the coefficients in the expansion of the same function with respect to a different set of Jacobi polynomials. They typically operate by computing the Chebyshev coefficients of an expansion and then applying a connection matrix or a series of connection matrices to obtain the coefficients in the desired expansion. The computation of the Chebyshev expansion can be carried out efficiently in O(n log n) operations via the nonuniform FFT (NUFFT) [12,13] and the application of each connection matrix can be performed in a number of operations which grows linearly or quasi-linearly via the fast multipole method (FMM) [3,14,15,5] or the fast eigendecomposition of semiseparable matrices [16,17,5].
The second class of algorithms, of which [18] is a major example, make use of butterfly methods and similar techniques for the application of oscillatory matrices. Most algorithms in this category, including that of in [18], require precomputation with O(n 2 ) running times. However, [19] describes an algorithm based on this approach whose total running time is O(n log 2 (n)), when both parameters a and b are in the interval (−1/2, 1/2). It makes use of the observation that the solutions of Jacobi's differential equation can be accurately represented via non-oscillatory phase and amplitude functions to apply the Jacobi transform rapidly. That certain differential equations admit non-oscillatory phase functions has long been known [20]; a general theorem was established in [21] and a numerical method for the computation of non-oscillatory phase functions was given in [22]. In the case of Jacobi's differential equation, there exists a smooth amplitude function M (a,b) (t, ν) and a smooth phase function ψ (a,b) (t, ν) such that whereP (a,b) ν andQ (a,b) ν are referred as the modified Jacobi functions of the first and second kind, respectively. They are defined in terms of the Jacobi functions of the first and second kinds P (see [20] for definitions) via the formulas where The normalization constant is chosen to ensure the L 2 (0, π) norm ofP j=0 is an orthonormal basis for L 2 (0, π). The change of variables x = cos(t) makes the singularities in phase and amplitude functions for Jacobi's differential equation more tractable. Obviously, there is no substantial difference in treating expansions of the type (1) and those of the form and we will consider the latter form here, more specially, the latter form with scaling by quadrature weights. Moreover, we will denote by {t k } k=1,...,n and {w k } k=1...,n the nodes and weights of the trigonometric Gauss-Jacobi quadrature rule obtained by applying the change of variables x = cos(t) to (2). Again, we refer to transforms which use the set of points {t k } k=1,...,n as uniform and those which use some other set points as nonuniform.
In [19], it was observed that the second summation in (8) could be interpreted as the application of a structured matrix that is the real part of a Hadamard product of a NUFFT matrix and a numerically low-rank matrix. This leads to a method for its computation which takes quasi-linear time [13,23,24]. Expansions in terms of the functions of the second kind {Q (a,b) j } ∞ j=0 can be handled similarly. This non-oscillatory representation method will be described in more detail in Section 2.2.
In this paper, we generalize the one-dimensional Jacobi polynomial transform of [19] to the multi-dimensional case. The transformation matrix in the multi-dimensional case is still the real part of a Hadamard product of a multi-dimensional NUFFT matrix and a numerically low-rank matrix in the form of tensor products. However, the numerical rank of the low-rank matrix might increase quickly in the dimension following the algorithm in [19] and hence the multi-dimensional NUFFT in [12,13] might not be sufficiently efficient. To obtain a viable multi-dimensional Jacobi polynomial transform, we reformulate the Hadamard product into a new one, which is the Hadamard product of an over-sampled Fourier transform matrix and a numerically low-rank matrix, and propose a fast randomized SVD to achieve near optimality in the low-rank approximation. This leads to an efficient implementation of two-dimensional and three-dimensional Jacobi polynomial transforms. Besides, as will be demonstrated numerically, the new method proposed in this paper is more robust than the method in [19] and works for a and b in the regime (−1, 1), other than the regime (− 1 2 , 1 2 ) in [19]. The remainder of the paper is organized as follows. In Section 2, we will first briefly introduce the fast SVD via randomized sampling and the theory of non-oscillatory phase functions. In Section 3, a variant of the fast one-dimensional Jacobi polynomial transform in [19] is given; it is followed by a description of our fast multi-dimensional Jacobi polynomial transform. In Section 4, numerical results which demonstrate the efficiency of our algorithm are described. Finally, we will conclude our discussion in Section 5.

Preliminaries
In this section, we will revisit the linear scaling randomized SVD introduced in [25] and the nonoscillatory phase and amplitude functions in [19] to make the presentation self-contained.

Approximate SVD via Randomized Sampling
For a numerically low-rank matrix Z ∈ C N ×N with an O(1) algorithm to evaluate an arbitrary entry, [25] introduced an O(N r 2 ) algorithm to construct a rank-r approximate SVD, Z ≈ U 0 Σ 0 V * 0 , from O(r) randomly selected rows and columns of Z, where * denotes the conjugate transpose, U 0 and V 0 ∈ C N ×r , Σ 0 ∈ C r×r is a diagonal matrix with approximate singular values in the diagonal part in a descent order.
Here, we adopt the standard notation for a submatrix in MATLAB: given a row index set I and a column index set J, Z I,J = Z(I, J) is the submatrix with entries from rows in I and columns in J; we also use ":" to denote the entire columns or rows of the matrix, i.e., Z I,: = Z(I, :) and Z :,J = Z(:, J). With these handy notations, we briefly introduce the randomized SVD as follows.
where P is the resulting permutation matrix and R = (r ij ) is an O(r) × n upper triangular matrix. Define the important column index set Π col to be the first r columns picked within the pivoted QR decomposition.
3. Randomly sample rq columns and denote their indices by S col . Let J = S col ∪ Π col . Perform a pivoted LQ decomposition of Z :,J to get where P is the resulting permutation matrix and L = (l ij ) is an m × O(r) lower triangular matrix. Define the important row index set Π row to be the first r rows picked within the pivoted LQ decomposition.
4. Repeat steps 2 and 3 a few times to ensure Π col and Π row sufficiently sample the important columns and rows of Z.

5.
Apply the pivoted QR factorization to Z :.Π col and let Q col be the matrix of the first r columns of the Q matrix. Similarly, apply the pivoted QR factorization to Z * Πrow,: and let Q row be the matrix of the first r columns of the Q matrix. 6. We seek a middle matrix M such that Z ≈ Q col M Q * row . To solve this problem efficiently, we approximately reduce it to a least-squares problem of a smaller size. Let S col and S row be the index sets of a few extra randomly sampled columns and rows. Let J = Π col ∪ S col and I = Π row ∪ S row . A simple least-squares solution to the problem gives M = (Q col ) † I,: Z I,J (Q * row ) † :,J , where (·) † stands for the pseudo-inverse.
Then the low-rank approximation of Z ≈ U 0 S 0 V * 0 is given by In our numerical implementation, iterating Steps 2 and 3 twice is empirically sufficient to achieve accurate low-rank approximations via Algorithm 2.1. Similar arguments as in [26] for a randomized CUR factorization can be applied to quantify the error and success probability rigorously for Algorithm 2.1. But at this point, we only focus on the application of Algorithm 2.1 to the fast Jacobi polynomial transform without theoretical analysis.
Note that our goal in the Jacobi polynomial transform is to construct a low-rank approximation of the low-rank matrix, i.e., Z ≈ U V * with U and V ∈ C N ×r , up to a fixed relative error ε, rather than a fixed rank. Algorithm 2.1 can also be embedded into an iterative process that gradually increases the rank parameter r to achieve the desired accuracy. We denote the rank parameter as r ǫ when it is large enough to obtain the accuracy ǫ.
When Z origins from the discretization of a smooth function Z(x, y) at the grid points {x i } 1≤i≤N and {y i } 1≤i≤N , i.e., Z ij = Z(x i , y j ), another standard method for constructing a low-rank factorization of Z ≈ U V * is Lagrange interpolation at Chebyshev grids in x or y. For example, let C µ := {µ i } 1≤i≤r denote the set of r Chebyshev grids in the domain of x, U ∈ R N ×r be the matrix consisting of the i-th Lagrange polynomial in x corresponding to µ i as its i-th column, and V ∈ C N ×r be the matrix such that V ij is equal toZ(µ j , y i ), where· denotes the conjugate operator, then Z ≈ U V * by the Lagrange interpolation. Usually, an oversampling parameter q is used via setting r = qr ǫ . Then a rank- The fast nonuniform FFT in [13] and the one-dimensional Jacobi polynomial transform in [19] adopted low-rank approximation via Lagrange interpolation to deal with the low-rank term in their Hadamard products of low-rank and (nonuniform) FFT matrices without an extra truncated SVD. In this paper, we propose to use the randomized SVD via random sampling to obtain nearly optimal rank in the low-rank approximation.

Non-oscillatory phase and amplitude functions
Given a pair of parameters a and b in − 1 2 , 1 2 , and a maximum degree of interest N max > 27, we revisit the fast algorithms in [19] for constructing non-oscillatory phase and amplitude functions ψ (a,b) (t, ν) and M (a,b) (t, ν) such that andQ The polynomials with t and ν out of these ranges can be evaluated by the well-known three-term recurrence relations or various asymptotic expansions; the lower bound of 27 was chosen to obtain optimal numerical performance.
Next, we present a few facts regarding the phase and amplitude functions related to Jacobi's differential equations. It is well known that the functionsP (a,b) ν andQ (a,b) ν satisfy the second order differential equation where We refer to Equation (16) as Jacobi differential equation. Following the derivation in [19], we can show that the pair {P } of real-valued, linearly independent solutions of (16) satisfies where W is the necessarily positive constant Wronskian of the pair {P with C an appropriately chosen constant related to σ 1 . Note that ∂ t ψ (a,b) (t, ν) > 0 since W > 0. Hence, the non-oscillatory amplitude function of {P Through straightforward computation, it can be verified that the square of the amplitude function satisfies the third order linear ordinary differential equation (ODE) Obviously, Hence, we can specify the values of M (a,b) (t, ν) 2 and its first two derivatives in t at a point on the interval [ 1 Nmax , π − 1 Nmax ] for any ν using various asymptotic expansion for P . Afterwards, we uniquely determine M (a,b) (t, ν) via solving the ODE (22) using a variant of the integral equation method of [12] (or any standard method for stiff problems).
To obtain the values of ψ (a,b) (t, ν), we first calculate the values of via (21). Next, we obtain the values of the functionψ defined viã To evaluate C, we first use a combination of asymptotic and series expansions to calculateP and C can be calculated in the obvious fashion.
The above discussion has introduced an algorithm to evaluate ψ (a,b) (t, ν) and To achieve a fast algorithm, we note that it is enough to conduct calculation on selected important grid points of with their values on the important grid points.
It is well known that a polynomial p(t) (here p(t) is eitherP can be evaluated in a numerically stable fashion using the barcyentric Chebyshev interpolation formula [27] where x 1 , . . . , x k are the nodes of the k-point Chebyshev grid on a sufficiently small interval (σ 1 , σ 2 ) (such that p(t) is a polynomial of degree at most k − 1) and Hence, it is sufficient to evaluate ψ (a,b) (t, ν) and M (a,b) (t, ν), which lead toP , at a tensor product of Chebyshev grid points in t and ν, and calculate them at an arbitrary location (t, ν) via a bivariate piecewise barcyentric Chebyshev interpolation.
More explicitly, let m ν be the least integer such that 3 mν ≥ N max , and let m t be equal to twice the least integer l such that π denote the 16-point piecewise Chebyshev gird on the intervals with M t = 15(m t − 1) + 1 points in total; let denote the nodes of the 24-point piecewise Chebyshev grid on the intervals with M t = 23(m ν − 1) + 1 points in total. The piecewise Chebyshev grids (30) and (32) form a tensor product Then ψ (a,b) (t, ν) and M (a,b) (t, ν) can be obtained via barcyentric Chebyshev interpolations in t and ν via their values at the tensor product restricted to the piece (α i , α i+1 ) × (β j , β j+1 ) that (t, ν) belongs to for some i and j.
Here we summarize the operation complexity of the whole algorithm above. A detailed discussion can be found in [19]. ODE solvers are applied for O (log (N max )) values of ν, and the solution of the solver is evaluated at O (log (N max )) values of t, making the total running time of the procedure just described O log 2 (N max ) .
Once the values of ψ

Fast Jacobi polynomial transforms
Without loss of generality, we assume that ψ (a,b) (t, ν), M (a,b) (t, ν), P can be evaluated in O(1) operations for any t ∈ (0, π) and ν ∈ [0, N max ]. A fast algorithm for rapidly computing Gauss-Jacobi quadrature nodes and weights in the Jacobi polynomial transform in (8) has also been introduced in [19]. We refer the reader to [19] for details and assume that quadrature nodes and weights are available in this section.

One-dimensional transform and its inverse
Here we propose a new variant of the one-dimensional Jacobi polynomial transform forP (a,b) j (t) in [19]. The new algorithm simplifies the discussion of the fast algorithm. The transform forQ is similar. For our purpose, we consider the n th order uniform forward Jacobi polynomial transform with a scaling, that is, calculate the vector of values  given the vector  of the coefficients in the expansion here, t 1 , . . . , t n , w 1 . . . , w n are the nodes and weights of the n-point trigonometric Gauss-Jacobi quadrature rule corresponding to the parameters a and b. The "nonuniform" forward Jacobi transform does not require trigonometric Gauss-Jacobi quadrature nodes and weights.
In the uniform transform, the properties of the trigonometric Gauss-Jacobi quadrature rule and the weighting by square roots in (35) ensure that the n × n matrix W J and J (a,b) n is the n × n matrix representing the sum in (37). Hence, the inverse transform is a simple matvec by the transpose of W J (a,b) n . In the nonuniform case, J (a,b) n is usually a very ill-conditioned matrix and the inverse transform requires solving a challenging system of linear equations, which will be discussed in a separate paper in the future.
It follows from (14) that the (j, k)-th entry of J since the part consisting of low-degree polynomials is computed via three-term recurrence relations or various asymptotic expansions. Let V Note that G (a,b) n is the real part of a discrete Fourier integral transform that can be evaluated with a quasi-linear complexity by the algorithm of [23] for applying Fourier integral transforms and a special case of that in [24] for general oscillatory integral transforms. Motivated by this observation, the algorithm in [19] constructs a low-rank matrix A (a,b) n ∈ C n×(n−28) whose (j, k)-th entry is Denote the n × (n − 28) nonuniform FFT matrix, whose (j, k)-th entry is as N n . Then where "⊗" denotes the Hadamard product. It can be immediately seen from (42) that J (a,b) n can be applied through a small number of NUFFTs [12,13] and a simple direct summation for V (a,b) n in a quasi-linear scaling of n.

More specifically, let a low-rank approximation of
where u j and v j are column vectors for each j. Then the transform from (36) to (35) can be approximated by the following sum, where α V is the subvector of α with the first 28 entries, α G is the subvector of α with the rest entries, D u denotes a diagonal matrix with a column vector u on its diagonal. The formula (44) can be carried out by r NUFFTs in O(rn log n) arithmetic operations. In practice, inspired by the NUFFT in [13], we can replace the Hadamard product of a low-rank matrix and a NUFFT matrix in (42) with a Hadamard product of a low-rank matrix and a DFT matrix, the matvec of which can be carried out more efficiently with a few numbers of FFTs.
Consider another matrix B (a,b) n ∈ C n×(n−28) whose (j, k)-th entry is defined via where [x] denotes the integer nearest to x. Then we have where F n is an n × (n − 28) matrix whose (j, k)-th entry is exp i(2π [t j n/2π] n (k + 27)) .
It's obvious that F n in (46) is a row permutation of an inverse DFT matrix. Note that the difference between the phase functions of A (a,b) n and B (a,b) n is less than π. Hence, B (a,b) n is also a low-rank matrix and we can apply the randomized SVD in Algorithm 2.1 to construct a low-rank approximation of B (a,b) n in O(r 2 n) operations. Suppose we have constructed an approximate rank-r SVD of B (a,b) n up to a desired accuracy using Algorithm 2.1; that is, suppose that we have computed the factorization where U ∈ C n×r , V ∈ C r×(n−28) , and S ∈ R r×r is a positive definite diagonal matrix. By rearranging the factors above, we have where u i and v T i denote the i th column vector of U S 1 2 and the i th row vector of S 1 2 V , respectively, and T denotes the matrix transpose. Once (49) is ready, we have where D u denotes a diagonal matrix with u on its diagonal. Formula (50) indicates that the Jacobi polynomial transform can be evaluated efficiently via r inverse FFTs, which requires O(rn log n) arithmetic operations and O(rn) memory. Compared to (44) used in the original fast Jacobi transform in [19], the number of inverse FFTs in (50) has been optimized. Thus Formula (50) would take fewer operations to compute a multidimensional transform. Note that, compared to the original method in [19], though our new method using (50) would cost more time to construct a low-rank approximation, the optimized rank of the low-rank approximation would accelerate the application of the multi-dimensional Jacobi transform. In addition, as we mentioned previously, the new method works in a larger range of a and b. In Section 4, we will provide numerical comparisons to demonstrate the superior of the new method over the method in [19].
The fast inverse Jacobi transform in the uniform case can be carried out in a similar manner, since W J (a,b) n is an orthonormal matrix. In fact, the inverse transform can be computed via where F T n is a permutation of the DFT matrix. Therefore, the inverse Jacobi polynomial transform can also be computed via r FFTs.
It is worth emphasizing that the fast algorithm introduced in this section also works for nonuniform Jacobi polynomial transforms since the low-rankness of B (a,b) n is independent of the samples in t and the quadrature weights in the uniform Jacobi polynomial transform. There will be numerical examples later to verify this. The inverse transform in the non-uniform case requires solving a highly ill-conditioned linear system, which will be reported in a separate paper in the future.

Two-dimensional transform and its inverse
In this section, we extend our algorithm to the two-dimensional case, using the new method developed in the last section. We only focus on the transform forP (a,b) ν (t); the one forQ (a,b) k−1 (t) is similar. We will adopt the MATLAB notation A(:) as a vector resulting from reshaping the matrix A into a vector. We also use similar notations as in Section 3.1, e.g., J n,· , F n,· , and W · denote corresponding matrices analogous to their counterparts in Section 3.1, respectively, with "·" specifying a variable x or y in the spatial domain. In the rest of this section, we always assume a low-rank approximation has been obtained by Algorithm 2.1 up to a desired accuracy for "·" as x or y.
Given locations {x i } i=1,···,n ⊂ (0, π) and {y i } i=1,···,n ⊂ (0, π), with no substantial difference, the forward and inverse two-dimensional Jacobi polynomial transforms arise in the following Jacobi expansion where α denotes an expansion coefficients matrix. The forward and inverse transform can be defined analogously as in the one-dimensional case. When both {x i } i=1,···,n and {y i } i=1,···,n are exactly the nodes of the trigonometric Gauss-Jacobi quadrature rule, the corresponding transform is referred as the uniform transform. With additional weight matrices W x and W y , the forward transformation n,y ) taking α(:) to (W x ⊙ W y )f (:) is orthogonal, where "⊙" denotes the Kronecker product. When {x i } i=1,···,n , {y i } i=1,···,n , and weights are not given by the trigonometric Gauss-Jacobi quadrature rule, we refer the corresponding transform as a non-uniform transform.
With the low-rank approximations in (52) available, we have V (a,b) n,y 0])α(:) where 0 is an n × (n − 28) zero matrix. In our numerical implement, since 0 does not contribute to the numerical result, there are just O(n) entries of α(:) needed to be applied in the first three terms above. Formula (54) indicates that 2D uniform Jacobi polynomial transforms can be evaluated via O((r x + r y )n) 1D inverse FFTs and O(r x r y ) 2D inverse FFTs, which results in O(r x r y n 2 log(n)) arithmetic operations and O(r x r y n 2 ) memories. We would like to emphasize that this algorithm also works for 2D non-uniform forward transform since it is actually a tensor form of the onedimensional transform. There will be numerical examples later to verify this.
The fast 2D inverse Jacobi transform in the uniform case can be carried out in a similar manner, since both W x J (a,b) n,x and W y J (a,b) n,y are orthonormal matrices. In fact, the 2D inverse transform can be computed via wheref (:) := (W x ⊙ W y )(W x ⊙ W y )f (:). In our numerical implement, since 0 does not contribute to the numerical result, we just need to compute O(n) entries of α(:) for the first three terms above. Therefore, 2D inverse Jacobi polynomial transforms can also be evaluated via O((r x + r y )n) 1D FFTs and O(r x r y ) 2D FFTs.

Three-dimensional transform and its inverse
In this section, we continue to extend our algorithm for 3D Jacobi polynomial transform and its inverse. Analogously, we just discuss the transform forP (a,b) ν (t) and the transform forQ k−1 (t) is similar. The notations in Section 3.2 will be inherited here.
Given locations {x i } i=1,···,n , {y i } i=1,···,n , and {z i } i=1,···,n in (0, π), with no substantial difference, the forward and inverse three-dimensional Jacobi polynomial transforms arise in the following Jacobi expansion where α denotes a three-dimensional tensor containing expansion coefficients. The forward and inverse transforms can be defined analogously to the three-dimensional case. When {x i } i=1,···,n , {y i } i=1,···,n , and {z i } i=1,···,n are all exactly the nodes of the trigonometric Gauss-Jacobi quadrature rule, the corresponding transform is referred as the uniform transform. Otherwise, we refer the corresponding transform as a non-uniform transform.
In the uniform transform, in order to take advantage of orthogonality, we consider the the tensor n,z ) that takes α(:) to (W x ⊙ W y ⊙ W z )f (:). Once the related low-rank approximations, have been obtained, we have n,y 0])(I n ⊙ F n,y ⊙ I n )(I n ⊙ D v y,i ⊙ I n )α(:)) n,z 0])(F n,x ⊙ F n,y ⊙ I n )(D v x,j ⊙ D v y,k ⊙ I n ))α(:) In numerical implement, since 0 does not contribute to numerical result, there are just O(n 2 ) entries of α(:) needed to be applied in the first seven terms above. Formula (59) indicates that 3D uniform Jacobi polynomial transform can be evaluated via O((r x + r y + r z )n 2 ) inverse FFTs, O((r x r y + r x r z + r y r z )n) 2D inverse FFTs and O(r x r y r z ) inverse 3D FFTs, which results in O(r x r y r z n 3 log(n)) arithmetic operations and O(r x r y r z n 3 ) memories. Again, we would like to emphasize that this algorithm also works for 3D non-uniform forward transform. There will be numerical examples later to verify this.
The fast 3D inverse (uniform) Jacobi transform in the uniform case can be carried out in a similar manner, since W x J n,z are all orthonormal matrices. In fact, the 3D inverse transform can be computed via wheref (: In our numerical implement, since 0 does not contribute to the summation above, we just need to compute O(n 2 ) entries of α(:) for the first seven terms above. Therefore, 3D inverse Jacobi polynomial transform can also be evaluated via O((r x +r y +r z )n 2 ) 1D FFTs, O((r x r y +r x r z +r y r z )n) 2D FFTs, and O(r x r y r z ) 3D FFTs.

Numerical results
This section presents several numerical examples to demonstrate the effectiveness of the algorithms proposed above. Section 4.1 provides a comparison of the original method in [19] named as method CHEB using (44) in Section 3.1, and our new method named as RS using (50) in Section 3.1 to demonstrate the superiority of our new method. In Section 4.2, we apply our new method to three examples corresponding to 1D, 2D, and 3D Jacobi polynomial transforms with parameters a = b = 0.25, respectively. All implementations are in MATLAB R on a server computer with 28    processors and 2.6 GHz CPU. We have made our code, including that for all of the experiments described here, available on GitLab at the following address: https://gitlab.com/FastOrthPolynomial/Jacobi.git In our numerical results, N denotes the number of grid points per dimension; "fac" and "app" stand for the factorization time for setting up the algorithm and the application time for applying the algorithm, respectively; "relerr" means the relative error of the fast matvec. The rank parameter r and the accuracy parameter ǫ in all low-rank factorization algorithms are set to ⌈2 log 2 (N )⌉ and 1e − 8, respectively, where ⌈·⌉ is the ceiling function. To make fair comparisons, the oversampling parameter q in Algorithm 2.1 for the RS method is set such that q × r is approximately equal to the number of piecewise Chebyshev grids in (30) used in the CHEB method. We perform the same experiment 10 times and summarize the average statistics in the following figures.

Comparison of CHEB and RS methods
4.1.1 Comparisons in the three-dimensional case for a and b in (− 1 2 , 1 2 ) First, we provide performance comparisons of the CHEB algorithm and the RS algorithm in the three-dimensional case with a = b = 0.25 to demonstrate that the RS is more efficient than the CHEB algorithm. In the one and two-dimensional cases, RS is also faster than CHEB but the speed-up is less obvious.
As we can see from the results summarized in Figure 1, the RS method provides a more compact matrix compression while keeping the compression accuracy competitive to that of the CHEB algorithm. The more compact compression results in faster set-up and application time for the Jacobi transform for all problem sizes.

4.1.2
Comparisons in the one-dimensional case for a and b in (−1, 1) \ (− 1 2 , 1 2 ) Second, we provide comparisons of the CHEB and RS algorithms when a and b are assigned as 0.6, 0.7, 0.8, and 0.9 in the range (−1, 1) \ (− 1 2 , 1 2 ), which is the range of a and b that has not been explored in [19]. Since the case when a and b are negative is similar to the case when a and b are positive, we only show the results for positive a and b. The numerical results are summarized in Figure 2 and they show that both CHEB and RS achieve desired accuracy when a and b are not close to 1; RS works but CHEB fails when a and b are close to 1. Hence, the newly proposed RS algorithm works in a larger range of a and b than the CHEB method in [19].

Performance of fast multi-dimensional transforms
Finally, we will present examples for two and three-dimensional Jacobi polynomial transforms using the RS algorithm. In either the two or three-dimensional case, there are numerical results for one uniform forward transform, one uniform inverse transform, and one nonuniform forward transform. a and b are all set to 0.25 in these examples. Figure 3 and 4 summarize the numerical results in 2D and 3D, respectively. In cases, no matter forward or inverse transform, uniform or nonuniform transform, 2D or 3D, the factorization and application time of our algorithm scale like O(N d log N ) or O(N d log 2 N ), where d is the dimension. The numerical rank gradually increases like O(log N ) as the problem size increases. The relative error of the fast algorithm is in line with the desired accuracy in the low-rank factorization.

Conclusion
This paper proposed a fast algorithm for multi-dimensional Jacobi polynomial transforms based on the observation that the solution of Jacobi's differential equation can be represented via nonoscillatory phase and amplitude functions. In particular, the transformation matrix corresponding to the discrete transform is a Hadamard product of a numerically low-rank matrix and a multidimensional discrete Fourier transform matrix. After constructing a low-rank approximation to the numerical low-rank matrix, the application of the Hadamard product can be carried out via O(1) fast Fourier transforms, resulting in a nearly optimal algorithm to compute the multi-dimensional Jacobi polynomial transform.
We proposed to apply the randomized SVD to construct low-rank factorizations, resulting in a faster Jacobi transform in high-dimensional spaces. Moreover, numerical experiments show that the new fast transform works for a larger class of Jacobi polynomials with parameter a and b in the interval (−1, 1) than the one-dimensional algorithm in [19].
For other values of a and b outside (−1, 1), the Jacobi transformation matrix is no longer purely oscillatory and hence the proposed method is not applicable. We will tackle this issue in the future via hierarchically dividing the transformation matrix into sub-matrices that are either purely nonoscillatory, which can be handled by fast low-rank factorization, or purely oscillatory, which can be applied via the fast algorithm in this paper.
The fast inverse in the nonuniform case is still an open problem. It involves a highly illconditioned linear system of equations. We are working on an efficient preconditioner for this linear system and will summarize our work in a separate paper.


extfieldsofstudy:
Computer Science
Mathematics

s2fieldsofstudy:
Computer Science
Mathematics

text:
Per-Flow Cardinality Estimation Based On Virtual LogLog Sketching

Flow cardinality estimation is the problem of estimating the number of distinct elements in a data flow, often with a stringent memory constraint. It has wide applications in network traffic measurement and in database systems. The virtual LogLog algorithm proposed recently by Xiao, Chen, Chen and Ling estimates the cardinalities of a large number of flows with a compact memory. The purpose of this thesis is to explore two new perspectives on the estimation process of this algorithm. Firstly, we propose and investigate a family of estimators that generalizes the original vHLL estimator and evaluate the performance of the vHLL estimator compared to other estimators in this family. Secondly, we propose an alternative solution to the estimation problem by deriving a maximum-likelihood estimator. Empirical evidence from both perspectives suggests the near-optimality of the vHLL estimator for per-flow estimation, analogous to the near-optimality of the HLL estimator for single-flow estimation.

today make it infeasible to spend much time on processing each element. According to [1] and [2], on a typical OC-768 2 backbone network link with 40 Gbps traffic speed, the time available to process each packet is at best about tens of nanoseconds, corresponding to no more than a hundred elementary operations.
Secondly, such a high traffic rate necessitates the need to use on-chip cache memory on network processors, in order to achieve real-time processing and maintain high throughput. But most on-chip caches on processors are made of SRAM, typically only a few megabytes [3]. A naïve approach that maintains a look-up dictionary to record all the distinct elements seen thus far is too costly, since it requires O(n) bits. In some applications, n can be on the order of billions.
Therefore, an ideal algorithm should be simple and quick in processing each flow element and be compact in memory usage. The good news is that it is usually not necessary to know the exact value of a flow's cardinality. Many approximation algorithms have been developed to explore the trade-off between estimation accuracy and space/timing efficiency.
Applications Examples: One important direction of application of cardinality estimation is network anomaly detection. Estan et al. [4] suggested three such applications: detecting port scans, detecting denial of service (DoS) attacks and estimating worm spread rate.
Let us take port scan detection as an example to illustrate the idea. A port scan is a probing process that sends service requests from a client to a range of port addresses of a server with the aim of finding an open port, which an attacker can take advantage of to find vulnerabilities on the server. A router can detect such a process by keeping track of packet flows and using cardinality estimation to measure the number of distinct port addresses attempted by each source. Any source that is trying to connect to an abnormally large number of port addresses in a short time interval should be suspected of conducting a port scan.
Other than network anomaly detection, we can apply flow cardinality estimation to many other problems that have a count-distinct nature by adapting the concept of a flow to contain other (maybe abstract) types of data. For example, Google might be interested in estimating the number of distinct users that query certain keywords in a period of time [3]. Here a per-keyword flow can be defined to be all the user IP addresses that query that keyword. Information about the cardinalities of such flows reflects the popularity of those keywords which might be helpful in the optimization of Google's database and search algorithm. Much research on cardinality estimation was motivated by database-related applications (e.g. [5] and [6]).

Motivation and contributions
State-of-the-art algorithms for cardinality estimation comprise two processes. First, the sketching process reads elements from the flow and stores useful information in a compact data structure. Second, in the estimation process, an estimator takes the recorded information as input and outputs the estimated cardinality of the flow.
The HyperLogLog (HLL) algorithm proposed by Flajolet et al. in [7] is near-optimal and has been widely adopted in many industry systems due to its simplicity and excellent performance; according to [3], it is the best existing algorithm for single-flow estimation. Roughly speaking, the HLL algorithm requires hundreds of bytes to make a fair estimation for a single flow with cardinality up to 4 × 10 9 . Details of the HLL algorithm are presented in Section 2.1.1.
However, in some applications the cardinalities of multiple flows need to be estimated at the same time -this is known as per-flow cardinality estimation. In this case, the rate of hundreds of bytes per flow is still too much in some important scenarios where the number of flows can be on the order of tens of millions and memory usage is critical. To this end, Xiao et al. [3] proposed the virtual LogLog algorithm based on memory sharing at register level, which can potentially bring down the memory cost from hundreds of bytes per flow to one bit per flow on average. Details of the sketching process of this algorithm are described in Section 2.2.1. The estimation process of this algorithm uses an estimator called vHLL, summarized in Section 3.2.1.
The main purpose of this thesis is to propose and investigate alternative estimators for the estimation process of the virtual LogLog algorithm in [3]: • The original vHLL estimator of [3] is based on the HLL estimator for single-flow estimation. We show that the HLL estimator can be generalized by a family of estimators parametrized by a value θ (Section 3.1); for HLL, θ = −1. The idea of this generalization can be easily applied to the vHLL estimator for per-flow estimation as well.
Although it is already known that for single-flow estimation, θ = −1 (i.e. HLL) is near-optimal [7], it is not clear whether this near-optimality at θ = −1 extends to the case of per-flow estimation. We provide empirical evidence to show that indeed θ = −1 (i.e. vHLL) is still optimal for per-flow estimation.
• We introduce an alternative approach to the estimation problem -a maximum-likelihood estimator. We find that this estimator does not significantly improve the estimation process, compared to the vHLL estimator.
Empirical evidence from both perspectives suggests the near-optimality of the vHLL estimator for per-flow estimation, on the basis of the virtual LogLog algorithm [3].

Thesis overview
The rest of this thesis is organized as follows. Chapter 2 provides the background knowledge and preliminaries to the estimators discussed in this thesis. Chapter 3 and Chapter 4 investigate alternative estimators of the virtual LogLog algorithm from the aforementioned two perspectives. Chapter 3 introduces the vLL θ estimator, a generalization of the vHLL estimator, and evaluates its performances for different values of θ. Chapter 4 derives the maximum-likelihood estimator and compares its performance with that of the vLL θ estimator. Chapter 5 summarizes the conclusions of the thesis and points to possible future works.
This chapter briefly overviews cardinality estimation algorithms and points to relevant references. The focus is on summarizing previous works that are essential to the understanding of the main body of this thesis: the LogLog [8] and HyperLogLog [7] algorithms for single-flow estimation and the virtual LogLog algorithm [3] for per-flow estimation.

Single-flow cardinality estimation
Given a single flow of elements x 1 , x 2 , x 3 , · · · , the goal of cardinality estimation is to estimate the number of distinct elements in the flow. There are two general approaches to this problem: sampling and streaming. A sampling-based algorithm collects a small sample of elements from the flow while bypassing other elements and infers the cardinality of the flow from sample information. A streaming-based algorithm, in contrast, processes every element in the flow. Most state-of-the-art algorithms are based on streaming. For a review of these two kinds of algorithms and a comparison between them, see [9].
A streaming-based algorithm for cardinality estimation has three integral components: • A hash function H, mapping each element x i of the flow into a hashed value H(x i ) of a chosen type. The purpose of H is generally twofold: filtering out duplicate elements (relying on the fact that multiple appearances of the same element have the same hash) and providing randomization.
• A compact data structure D, capturing certain low-dimensional statistical information of the flow's cardinality by recording features of the hashed values. In some literature, the data structure is known as a sketch, meaning it is a concise summary of the flow. The process of updating the sketches is called sketching.
• An estimator E that takes the content of D as input and outputs the estimated cardinality, thus n = E(D).
Usually we want the estimator to be unbiased (i.e. E[ n|n] = n). One commonly used metric for the evaluation of the performance of an estimator is the relative standard error √ Var( n) n . For two unbiased estimators using the same amount of memory, the one with a smaller relative standard error is better. For a detailed classification and comparison of existing single-flow cardinality estimation algorithms, see [10,11]. We take a closer look at the LogLog and HyperLogLog algorithms in Section 2.1.1.

LogLog and HyperLogLog
The HyperLogLog (HLL) [7] algorithm is the best-known algorithm in practice, due to its simplicity and good performance. HLL is a successor of LogLog [8]: the two algorithms use the same hash function and data structure (therefore have the same sketching process) and only differ in their choice of estimators.
The data structure of the LogLog/HLL algorithm is a single array of k registers (counters), denoted as S, where S[i] refers to the i th register in the array (or the value stored in it, depending on the context). Suppose we have a hash function H that maps each element of the flow to a sufficiently long bit string. Define a function ρ that takes a bit string and outputs the position of the leftmost 1-bit of the string, e.g. ρ(000010 · · · ) = 5. The sketching process of LogLog/HLL is summarized as follows.
LogLog/HLL sketching process: 1. Initialize all the registers in S with value 0. Let l = log 2 k.
2. For each element x in flow: . Hash x into a sufficiently long bit string.
(ii) Divide the bit string into two parts: Use the prefix l bits to select a register in S.
Use the remaining bits to update the selected register. Here are a few points to note about the above sketching process. First, k is chosen to be some integer power of 2 so that l is an integer and j ∈ {0, 1, · · · , k − 1}. Assume that for a random element x, H(x) is a random bit string with independent and uniform bits. Then j can be regarded as a uniform random variable in {0, · · · , k − 1}, which means that the k registers in S are equally likely to be selected for updating for a randomly given element. Second, the values of the registers in S after the sketching process should not be affected by duplicate elements, nor by the order of the elements appearing in the flow. This is guaranteed by hashing and the max operation in step (iii).
Third, how large can a register's value get? With the previous assumption on the independence and uniformity of the bits in the hash output string, for a random element x we have ρ(q) = i with probability 1 2 i (i ≥ 1). That is, ρ(q) can be regarded as a geometric random variable with parameter 1/2. By the same assumption, we expect about n k distinct elements distributed to any particular register S[j]. Suppose that this number n k is exact, then at the end of the sketching process, S[j] is just the maximum of n k independent Geo(1/2) random variables. According to [8], previous study has shown that the expectation of S[j] is close to log n k with a small additive bias. Therefore, on average each register needs about log log n k = O(log log n) bits to store its value -this is why such a register is also known as a LogLog sketch. Based on this discussion, a rough estimator for n can be n = k2 S[j] . Following this idea, the LogLog estimator [8] replaces S[j] with the average of the k register values: n = LL(S) ηk2 where η is a suitable bias correction factor. The HLL algorithm [7] uses a different estimator: with a different bias correction factor γ. Note that 2 is commonly known as the geometric mean of the values 2 We will introduce two other estimators for single-flow estimation in Section 3.1 (LL θ estimator) and Section 4.4 (maximum-likelihood estimator), respectively.
It has been shown in [8] and [7] that the LogLog and HLL estimators are asymptotically approximately unbiased in the sense that, as n → ∞, E[ n] n is very close to 1 with a practically negligible fluctuation. It has been shown that as n → ∞, the bias correction constants are independent of n. They do depend on k, however; but as k gets large (e.g. k ≥ 64), they can be safely replaced by constants for all practical purposes. Practical values of η and γ are 0.39701 and 0.7213 respectively.
The relative standard error √ Var( n) n is approximately 1.30 √ k for the LogLog estimator [8] and approximately 1.04 √ k for the HLL estimator [7], as n → ∞. Therefore, using the same number of registers (same amount of memory), HLL's estimate is more accurate than that of LogLog. With the HLL algorithm, we can achieve an estimate accuracy (in terms of relative standard error) of about 5% by choosing k = 512. If we use 5 bits for each register (for measuring cardinalities up to 2 2 5 ≈ 4 × 10 9 ), then in total we need 5k = 2560 bits = 320 bytes for one estimation.
According to Section 4 of [7], the HLL algorithm is near-optimal, in the sense that its relative standard error 1.04 √ k is quite close to the lower bound 1 √ k for a wider class of algorithms based on order statistics. For a more detailed discussion of this lower bound, see [2].
One problem of the HLL algorithm is that it is highly biased and inaccurate for the estimation of small cardinalities. Methods to remedy this flaw are considered in Section 4 of [7] and further in [6].

Per-flow cardinality estimation with memory sharing
In many real-world applications we need to estimate the cardinalities of multiple flows at the same time. For example, consider a stream of data packets from many flows observed by a network monitor device (such as a router). Let each packet be abstracted as a 2-tuple (f, x), where f is the IP source address of the packet (the flow ID) and x is the destination address (i.e. an element of the flow). The goal is, at the end of the measuring period, to estimate the number of distinct destination addresses from each given source address, i.e. the cardinality of each flow. For example, if the stream of packets is (A, 2), (C, 9) An immediate idea to solve this problem is to allocate a separate block of memory for each flow and use any of the existing single-flow cardinality estimation algorithms for each flow. Since we do not know the cardinalities of the flows beforehand, it is inevitable to allocate the maximum amount of memory for each flow; with the HLL algorithm, we still need hundreds of bytes per flow. However, in many applications most of the flows have small cardinalities.  In view of this waste of memory, some algorithms have been proposed that allow memory to be shared among flows. Such an algorithm usually combines the following components/ideas: • A memory pool, which is the actual source of memory for all estimations.
• A virtual data structure for each flow. This data structure does not physically exist but is logically constructed by (often randomly) pointing to memory units in the pool. Consequently, difference flows may share some common parts of the memory.
• In the sketching process, for each incoming element, the algorithm applies the sketching process of an existing single-flow estimation algorithm to update the virtual data structure allocated for the corresponding flow.
• After the sketching process, the algorithm uses an estimator to estimate the cardinality of any given flow (with its ID). The estimator needs to take into account the interference among flows due to memory sharing. This estimation process can be done off-line (in comparison to the online sketching process).
Such algorithms may differ from each other in any of the above aspects. A review and comparison of some per-flow cardinality estimation algorithms based on memory sharing can be found in [3]. All these algorithms, however, share memory at bit level (meaning that the basic unit of the memory pool and of the virtual data structure is a bit), which [3] claims to be inherently too noisy.
Another algorithm that shares memory at bit level, which [3] did not mention, is the virtual FM sketches algorithm proposed in [13]. The algorithm constructs multiple virtual bit arrays from a bit pool for each flow and relies on the sketching process of the FM sketches 2 algorithm for single-flow estimation. The algorithm adopts a maximum-likelihood estimator in the estimation process based on bit patterns in the sketches. We will also consider a maximum-likelihood estimator solution in Chapter 4, but for a different sketching scheme. The performance reported in [13] shows that with memory cost of 1 bit per flow, the algorithm can estimate flows with cardinalities up to 3000 and average flow cardinality about 2.5; it can achieve relative standard error of about 20% for flows with cardinalities less than 500 and 10% for flows with cardinalities in the range 500 to 3000. For relatively small flows with cardinalities in the range 100 to 500, this algorithm may outperform the virtual LogLog algorithm with the vHLL estimator [3]; but it does not appear to be as competitive for larger ranges of flow cardinality.
In view of the inherent drawback of bit-level sharing, [3] proposed the virtual LogLog algorithm based on memory sharing at register level and showed through experiments that such an algorithm outperforms previous ones. The data structure and sketching process of this algorithm are presented in the following subsection.

Virtual LogLog register sharing and sketching process
Virtual Data Structure: The virtual LogLog algorithm [3] keeps a memory pool in the form of a register array, denote as R. Denote the size of the array by m, which is usually a very large number. R[j] refers to the j th register in the array (or the value stored in the register, depending on the context). For each flow f , we form a virtual data structure (virtual register array) denoted as R f , which is a logically constructed array of k registers with the i th register denoted as R f [i]. The registers of R f are randomly selected from R by using k independent hash functions G 0 , G 1 , · · · , G k−1 , each mapping the flow ID uniformly to an integer in {0, · · · , m − 1}, i.e.
The k hash functions can be implemented by a single master hash function G as follows: where "|" is the concatenation operator. It should be emphasized that R f does not need to be physically constructed (thus it is "virtual"). A simple example is shown in Figure 2.2 to illustrate the concepts. In the example, say we want to update the third register of flow f 2 , i.e. R f 2 [3], what actually will be updated is R[5] -we do not even need to know where the other registers in R f 2 are. There are three flows and three corresponding virtual register arrays: A caveat should be pointed out here, which the original paper [3] omitted: it could happen that two logically distinct registers of a virtual register array are mapped from the same physical register in the pool, i.e. G i (f ) = G j (f ) for some i = j. Certainly we wish to avoid this situation, but it is tolerable if the number of such "collisions" is very small.
To see how likely it is for such a "collision" to occur, we can consider a bins-and-balls analogous problem: suppose we have k balls and n bins and we throw the balls sequentially, independently and uniformly at random into the bins. In the end how many bins will contain more than one balls? If k is reasonably large, then the number of balls distributed at each bin can be approximated by an independent Poisson random variable X with mean k m (see Chapter 5.4 of [14]). So we have where the last approximation assumes k m. Therefore out of the m bins, we expect to have approximately k 2 m bins that hold more than one ball. Ideally we want k 2 m to be small, which means with high probability the registers in a virtual register array are all mapped from distinct physical registers in the pool. This factor should be included in the design consideration.
Sketching Process: The sketching process of virtual LogLog is almost identical to that of LogLog/HLL (recall from Section 2.1.1): for each packet (f, x) in the stream, we process x to obtain j (the register selector) and ρ(q) (to be compared with the selected register's value); except that the last step (iii) becomes That is, we are treating the virtual register array as the actual register array. By combining (2.3) and (2.4), the above expression can be re-written as Again, it shows that updates are actually made in the physical registers in R. It should now be clear why the algorithm is called virtual LogLog: it is based on virtual register arrays, where each register is a LogLog sketch. After the sketching process, we obtain m register values. The remaining problem is to infer flow cardinalities from these register values: this is the estimation process. We will consider two kinds of estimators in Chapters 3 and 4, respectively.

A per-flow estimator performance metric -Weighted square error
Before we discuss specific estimators for per-flow cardinality estimation, we consider a metric based on weighted square error to evaluate the performance of any given estimator.
Suppose that the incoming data stream contains M flows, with cardinalities n 1 , n 2 , · · · , n M . Given any estimator E, suppose its corresponding estimates for the flows' cardinalities are n 1 , n 2 , · · · , n M . The weighted square error of estimator E is defined as where w(·) is a weight function mapping the cardinality of a flow to a positive number. For two estimators using the same amount of memory, the one with a lower WSE is better.

Weight function
The choice of the weight function depends on the specific application. For example, if each flow is considered equally important regardless of its cardinality, then we can simply let w(n) = 1 for all n. In many other situations, large flows are considered to be more important than small flows, then we want w(n) be an increasing function in n. The weight function used in this thesis is presented and explained as follows. First, assume that the cardinality of a randomly chosen flow can be modeled as a random variable N with pmf p N (n) Pr(N = n). We use the following weight function: . (2.7) The motivation for using this weight function is explained here. The integral is an approximation of the total weight of flows whose cardinalities are in the range [m 1 , m 2 ). With the weight function in (2.7), we have w(n)p N (n) = 1 n ; it is easy to verify that in this case the integral for > 0 is independent of n * . If we plot the curve w(n)p N (n) as a function of n on log scale of n, then the area under the curve should be approximately the same in each decade interval: [1,10), [10,100), [100, 1000), [1000, 10000), etc. In another word, by choosing such a weight function, we put approximately the same total weight to the aggregate of flow cardinalities in each of these intervals.

Zipf model
The weight function in (2.7) can be applied to any given flow cardinality distribution. The particular distribution used for simulation in this thesis is the Zipf distribution. A random variable N following the Zipf(π, n max ) distribution has the following pmf: where n max is an upper bound of the flow cardinality, π is a parameter that controls the shape of the Zipf distribution (π > 0) and C is a constant such that nmax n=1 p N (n) = 1. The adoption of the Zipf model is motivated by the fact that Zipf's law underlies many Internet applications (see [15]). In Figure 2  With N ∼ Zipf(π, n max ), by (2.7) and (2.8), we have where the constant C can be omitted (set to 1) because when we compare the weighted square errors of two estimators, we only care about their relative values, which are unaffected by a constant factor. We remark that the specific estimators we will discuss in Chapters 3 and 4 do not reply on any particular distribution model. The Zipf model here is for two purposes: to complete the definition of the weighted square error for the performance evaluation of estimators and to generate random flow cardinalities for simulation; both purposes are independent of the estimation process.
Finally, we make a note on experimental evaluation of per-flow estimators. Given the close match between the real trace data and simulated trace data with a Zipf distribution, we will evaluate the estimators using simulated trace data. We generated 100 simulated trace files. Each trace file contains 10 6 flows, with the cardinality of each flow randomly and independently generated according to a Zipf(2.25, 10 5 ) distribution. A flow is presented in the file as a collection of packets with the same source address and distinct destination addresses; different flows have different source addresses (flow IDs). For one experiment, we process one such trace file and estimate the cardinalities of all the 10 6 flows in it. For a given estimator, we perform 100 independent experiments and evaluate its performance based on the statistics averaged over the 100 experiments. All experiments are performed with m = 200, 000 and k = 512. If each register uses 5 bits, this setting uses 10 6 bits to measure the cardinalities of 10 6 flows, leading to one bit per flow on average.

Chapter 3 VIRTUAL LOGLOG ESTIMATOR WITH PARAMETER θ
The sketching process of the virtual LogLog algorithm has been described in Section 2.2.1.
After the sketching process, we can offload the register values from the network measuring device for off-line query. The following per-flow estimation problem is to be solved: In this chapter, we first introduce LL θ , a family of generalized LogLog estimators parameterized by θ, for single-flow estimation. Then, with the similar idea from LL θ , we propose vLL θ , a family of generalized virtual LogLog estimators parameterized by θ, for per-flow estimation.

The LL θ estimator for single-flow estimation
The LL θ estimator is a class of estimators parameterized by θ, which unifies and generalizes the LogLog and HLL estimators described in Section 2.1.1.

Generalized mean
We start with the concept of generalized mean. Given k positive numbers x 1 , x 2 , · · · , x k , their generalized mean parameterized by θ is defined for nonzero θ by In the case of θ = 0, we let which is in fact the limit of A θ (x 1 , · · · , x k ) as θ → 0. Note that

Unification and generalization of LogLog and HLL estimators
Recall from Section 2.1.1, the LogLog and HLL estimators use the geometric mean and harmonic mean, respectively. Based on the generalized mean notation, we attempt to unify these two estimators by proposing the LL θ estimator: where ξ is a suitable coefficient. It is desirable if there exists a value of ξ to make the estimator approximately unbiased and for which we can identify such value. Specifically, ξ should not depend on n. We already know that for θ = 0 (i.e. LogLog [8]) and θ = −1 (i.e. HLL [7]) such values of ξ exist (by letting ξ = η and ξ = γ respectively). But we do not know if it is true for other values of θ. The analysis of this estimator for general θ is difficult. 1 We resort to simulations to explore empirical evidence of the existence of this coefficient ξ.
, · · · , 2 S[k−1] for short. In   From the plots we see: • For the same value of θ and k, n kA θ is almost a constant when n gets large (say n > 10000).
• For the same value of θ and n, n kA θ is almost a constant for large k (256, 512 or 1024). We conclude that, at least for θ ∈ [−3.0, 1.0], there exists a value for the aforementioned coefficient ξ that approximately depends only on θ for large n and k. To reflect this dependence of ξ on θ, we denote this coefficient as ξ θ instead.
To investigate the relationship between ξ θ and θ,  We observe an approximately linear relationship between ξ θ and θ by which is plotted as the red line in Figure 3.2. The LL θ estimator is now completely defined by replacing ξ in (3.4) with ξ θ specified in (3.5). The claim that this family of estimators unifies the LogLog and HLL estimators can be verified by checking that ξ 0 ≈ η (for LogLog) and ξ −1 ≈ γ for (HLL). A natural question one may ask next is if there exists an optimal value of θ for the LL θ estimator. To answer this question, we empirically calculate the relative standard error of the LL θ estimator for different values of θ (and for selected k's). The results are plotted in  According to the experiment results, for each fixed value of k, the optimal value of θ for the LL θ estimator is either −1 or −0.9; in the cases where LL −0.9 is the optimal, the difference between LL −1 and LL −0.9 , in terms of relative standard error, is negligible. Recall the LL −1 estimator is the same as the HLL estimator. This verifies the claim that the HLL is near-optimal [7], which we briefly described in Section 2.1.1.

The vLL θ estimator for per-flow estimation
With the LL θ estimator defined, we now introduce the vLL θ estimator for per-flow estimation. We start with a high-level idea of the estimator and the motivation for introducing the parameter θ here.
Just like for single-flow estimation where we infer the cardinality of a flow from its register array, for per-flow estimation we can infer the cardinality of a given flow f from its virtual register array R f . We can directly apply the LL θ estimator on R f to give an estimate of the total number of distinct elements distributed to R f . But this estimate involves the noise brought by other flows that have registers shared with flow f ; so we probably need a different bias correction coefficient to calibrate this rough estimate. While for the single-flow case θ = −1 is near-optimal, it is not immediately clear whether θ = −1 with a suitable bias correction coefficient is still near-optimal in the per-flow case.
Recall that A θ with a lesser value of θ is more robust to abnormally high values in obtaining the mean. For per-flow estimation, due to register sharing, large flows can cause much noise to some small flows by causing abnormally high register values. Therefore, we speculate that an estimator with a lesser value of θ might work better in the per-flow case due to its robustness against bursty noise from large flows. This motivates our investigation of the vLL θ estimator, introduced in the rest of this section.

Review of vHLL estimator and the generalization to vLL θ estimator
We start with a review of the virtual HyperLogLog (vHLL) estimator introduced in [3]. Let n T be the total aggregate cardinality of all flows, i.e. the sum of the cardinalities of all flows in the packet stream. Let n f + be the total number of distinct elements distributed to the register array R f , which include the distinct elements from flow f and those from other flows (we call noise elements). Suggested in [3] is the following approximate relationship between n f , n f + and n T : which can be interpreted as: the number of noise elements received by R f is the total number of elements from flows other than f scaled by k m -the ratio of the number of registers in R f to the number of registers in R. The assumption here is that noise elements are roughly uniformly distributed to all the m registers in the pool, which is a fair approximation when the number of flows and the number of registers for each flow are both sufficiently large [3]. Rearranging (3.6) gives us Now the values of n T and n f + are not directly available so we need their estimates n T and n f + . It has been shown in Section 6.1 of [3] that if n f + and n T are close to n f + k m (n T − n f ) (which approximately equals to n f + by (3.6)) and n T respectively, then n f is an approximately unbiased estimator of n f .
There are two possible methods to estimate n T . The first method is to treat n T as the cardinality of a grand flow -the flow containing all the distinct packets in the stream. In another word, if a packet is abstracted as a (f, x) pair, to estimate n T is to estimate the total number of distinct f − x pairs (e.g. source-destination pair) in the stream. This is a singleflow cardinality estimation problem and, as discussed in Section 2.1, can be solved using the HLL algorithm with an additional few hundreds of bytes, which is negligible compared to the memory for main register pool R. The second method is to use n T = HLL(R) as a rough estimator, based on the assumption that all the elements of the grand flow are distributed approximately uniformly over R. Either method works fine in practice. So n T is relatively easy to obtain. Now, for n f + , the vHLL estimator applies the HLL single-flow estimator on R f to obtain an estimate, i.e. (3.8) Combining the ideas above, the vHLL estimator for flow f is summarized as: We consider a generalization of the vHLL estimator by replacing the above harmonic mean A −1 with a generalized mean A θ . More specifically, given any value θ, let be a rough estimate of n f . Then we calibrate this rough estimate by suitable additive and multiplicative constants α θ and β θ to obtain a better estimate, i.e.
We discuss how to set the values of α θ and β θ in Section 3.2.2.

How to set α θ and β θ
For a given θ, one way to set the values of α θ and β θ is through empirical error minimization, explained as follows. Recall from the end of Section 2.3.2, we have 100 simulated trace files, each can be considered as a training sample. Suppose that in one training sample we have M flows with cardinalities n 1 , n 2 , · · · , n M , and by the vLL θ estimator specified in (3.11) we obtain corresponding estimates n 1 (θ), n 2 (θ), · · · , n M (θ). Then we can find values of α θ and β θ that minimize the weighted square error defined in Section 2.3 on this training sample. That is (n i − (α + β n i (θ))) 2 w(n i ), (3.12) where we let w(n i ) = n π−1 i and π = 2.25. The solution to the above minimization problem is standard: .
Since we have 100 training samples, we can obtain the values of α θ and β θ for each of the samples using the above method and use their average values for the estimator. The values of α θ and β θ obtained in this way are optimal in the sense that they minimize the weighted square error of the training samples (i.e. the empirical error ). The drawback of this approach is obvious. It only works well for the training samples generated according to a specific Zipf distribution; there is no performance guarantee for other data sets. Also, it relies on the specific definition of the weighted square error (including the weight function and the cardinality distribution model), which may not be universal for all applications. We suggest the following alternative and practical choice of α θ and β θ for the vLL θ estimator, which does not reply on any assumption on weight function or cardinality distribution model: where ξ θ can be calculated using (3.5) and n T can be calculated using either of the two methods mentioned in Section 3.2.1 (e.g. n T = HLL(R)). Note that this is a direct generalization of the vHLL estimator in (3.9). Also note that in this case α θ does not depend on θ.
In Figure 3.4, we plot the empirical values of α θ and β θ obtained by these two different methods through experiments on the 100 training samples. The plots show that our suggested values of α θ and β θ are good as they are close to the optimal values of α θ and β θ for these training samples.

Summary of the vLL θ estimator
We summarize the vLL θ estimator here. For a given θ ∈ [−3.0, 1.0], the vLL θ estimator estimates the cardinality of a given flow f by and n T is an estimate of the total aggregate cardinality of all the flows and can be obtained by calculating HLL(R) in practice. Note that it could happen that Equation (3.15) produces an estimate n f < 1, which is impossible in practice. In such case, we can simply reset n f = 1.

Experimental performance evaluation
We ran experiments on the 100 simulated trace files to evaluate the performance of the vLL θ estimator for different values of θ, in terms of the weighted square error defined in Section 2.3. The results are plotted in Figure 3.5. Results for θ ≥ 0.5 are not plotted in the graph because the values are comparatively too large (which means the corresponding vLL θ estimators are bad). According to the experiment results, the vLL −1 estimator has the best performance, which refutes our speculation at the beginning of Section 3.2 that a value of θ lesser than −1 might be optimal for the per-flow estimation case.
A possible explanation of this phenomenon is as follows. With the weight function defined in Section 2.3 and Zipf(2.25, 10 5 ) used for flow cardinality distribution, we have that w(n) ∝ n 1.15 . This means we put much larger weights on large flows compared to small flows. Therefore the weighted square error of an estimator is largely determined by how well the estimator estimates large flows. But the influence of bursty noise on large flows is much weaker than that on small flows; in another word, the signal-to-noise ratio for the cardinality estimation of a large flow is larger than that for small flows. Therefore, as we have discussed for the single-flow case where θ = −1 is near-optimal, for per-flow case θ = −1 is also likely to give good estimates for large flows and hence results in a smaller weighted square error.

Chapter 4 MAXIMUM LIKELIHOOD ESTIMATOR
In this chapter we derive an alternative approach to the same per-flow estimation problem stated at the beginning of Chapter 3: a maximum-likelihood estimator.

Formulation
We model the values of the k registers in R f after the sketching process as a random vector where Z f,i ∈ {0, · · · , r max }. The r max is the maximum value that can be stored in a register. For example, if the register has 5 bits, then r max = 2 5 − 1 = 31. Register values R f = (R f [0], · · · , R f [k − 1]) is an instance of the random vector Z f . Let p n (z 0 , · · · , z k−1 ) be the pmf of Z f (i.e. the joint pmf of the k random variables) given that n f = n, i.e. it is the likelihood function of n: p n (z 0 , · · · , z k−1 ) Pr Z f = (z 0 , · · · , z k−1 ) n f = n . (4.1) The maximum likelihood estimator finds the value of n that maximizes the likelihood function based on the instance of register values in R f . An equivalent formulation is to maximize the natural log of the likelihood function The maximum-likelihood estimator of flow f 's cardinality is then given by n f,ML = arg max n L f (n), n ∈ {1, 2, · · · , n max }, where n max is an upper bound of the cardinality of a flow.

Derivation of the likelihood function
In this section we derive an analytical expression for the likelihood function p n (z 0 , · · · , z k−1 ).

Outline of derivation
Suppose we have two copies of the register pool R: R (1) and R (2) , both initialized with 0's. Consider the following process: 1. Repeat the virtual LogLog sketching process described in Section 2.2.1 for all elements, except those from flow f , with register pool R (1) .
2. Repeat the sketching process only for elements from flow f with register pool R (2) .
3. Merge the register values in R (1) and R (2) by a register-wise max operation.
The register values in the merged pool as described above are exactly the same as those in R after the original sketching process. This is because the register values at the end of the sketching process are not affected by the order of arrival of the elements in the stream. Let R (1) f ) denote the virtual register array for flow f constructed by selecting registers from R (1) (R (2) ), whose indices in R (1) (R (2) ) are the same as those of R f in R. Then define the following random variables: Z the value of an arbitrary register in R (1) .
W n the value of an arbitrary register in R (2) f , conditioned on n f = n.
Z n the value of an arbitrary register in R f , conditioned on n f = n.
Z represents the background noise in our estimation caused by elements from flows other than f . W n represents the impact of elements from flow f itself. We model the values of the k corresponding registers in R (1) f and R (2) f also as random vectors, respectively, i.e.
Then we have f,1 , · · · , Z (1) Here is an outline of the derivation of an analytical expression for p n (z 0 , · · · , z k−1 ): 1. Find the distribution of Z and argue that Z (1) f,1 , · · · , Z (1) f,k−1 are independent. We will see that it is difficult to obtain the exact distribution of Z , but there exists a convenient and close approximation to it.
2. Find the distribution of W n and argue that Z f,k−1 are independent under the assumption that n is large.
3. With the distribution of Z and W n , find the distribution of Z n by (4.4).

Distribution of Z and independence of Z
(1) f,0 , · · · , Z (1) f,k−1 One possible approach to obtain an approximate expression of the distribution of Z is as follows. First assume that the total number of flows (excluding flow f ) is known and a flow's cardinality can be modeled as a random variable (e.g. Zipf). Then we can model the number of distinct elements distributed to an arbitrary register by a random variable Y and calculate the distribution of Y by adding up the influences from all the individual independent flows on that register. This calculation involves a high-dimensional convolution, which is computationally costly. We can perhaps avoid the convolution by appealing to the central limit theory for approximation. However this approximation gives us a continuous distribution for Y , from which we need to derive the CDF of Z that only takes discrete values (and really concentrates on only a few values as we shall see soon). This calculation may be highly inaccurate and is somewhat complicated. Moreover, in practice we usually do not know the total number of flows beforehand and it is undesirable that our estimator relies on any particular cardinality distribution model.
Another approach to approximate the distribution of Z , which circumvents the above difficulties and complications, is the following. Define a random variable Z: Z the value of an arbitrary register in R.
We claim that, for any practical purpose, the distribution of Z can be directly approximated by the distribution of Z, i.e.
where F Z and F Z are the CDFs of Z and Z, respectively. This approximation is based on the following two considerations: • Register values in R (1) differ from those in R at no more than k registers, because flow f can only affect k registers. If k m, the difference in the CDF of Z compared to that of Z is small.
• Assume there are many flows in the stream and n f is small compared to the total cardinalities of all other flows. As packets are being distributed to a register, it becomes harder and harder for the register's value to further increase, because it requires an much less likely (with geometrically decaying probability) hashed value to occur.
The analytical distribution of Z is also difficult to find. However, we can obtain an empirical distribution of Z directly from the register values in the pool at the end of the sketching process: In Figure 4.1 we show a sample distribution of Z generated using real trace data (the same data used for plotting Figure 2.1). From the figure we can see that the values of most of the registers are centered around 2 to 7.
Since Z is the distribution of an arbitrary register value in R (1) , including the k registers in R (2) f . We have Z (1) f,k−1 ∼ Z . Since the k registers are randomly selected from the pool with replacement, their independence is guaranteed by a good choice of the hash function G.

Distribution of W n and independence of Z
(2) f,0 , · · · , Z (2) f,k−1 To find the distribution of W n , we first define a random variable X n : X n the number of distinct elements distributed to an arbitrary register in R (2) f , conditioned on n f = n.
It is easy to see X n ∼ Binom n, 1 k . The number of distinct elements distributed to each of the k registers in R (2) f has the same distribution as X n , but they are not independent of each other (because they sum to n). However, if n is large, we can approximate n by a Poisson random variable with mean n. Under this Poissonization approximation, X n ∼ Poisson n k and the number of distinct elements distributed to each of the k registers in R (2) f will then be independent (see Chapter 5.4 of [14] for this Poissonization trick).
As discussed in Section 2.1.1, W n can be regarded as the maximum of X n i.i.d. geometric random variables with parameter 1/2. Given that X n ∼ Poisson n k , the CDF of W n can be calculated as follows. For 0 ≤ i ≤ r max , Under the Poissonization approximation, Z f,0 , Z f,1 , · · · , Z (2) f,k−1 are independent and have the same distribution with W n . We emphasize that this independence follows from the assumption that n is large.

Distribution of Z n and independence of
Recall that, by construction, Z n = max{Z , W n }. Therefore the CDF of Z n is: where F Z (i) and F Wn (i) are the CDF of Z and W n respectively. We have used that Z and W n are independent of each other and the distribution of Z can be approximated by that of Z. Hence, the pmf of Z n can be approximated by where with Poissonization of n. (4.9) We have argued the independence of Z (1) f,k−1 and the independence of Z (2) f,i by construction. So Z f,0 , · · · , Z f,k−1 are also independent of each other under the same conditions. The consequence of this independence is that we can factor the k-variable likelihood function p n (z 0 , · · · , z k−1 ) into the product of k single-variable pmfs, based on the distribution of Z n . That is, under this approximation, the log likelihood function in (4.2) can be re-written as (4.10)

Implementation
It would be nice if we could obtain a closed-form analytical expression of the value of n that maximizes the log-likelihood function L f (n). Unfortunately it turns out to be difficult. In order to search for this value, we explore properties of L f (n) that might be helpful.

Concavity of the log-likelihood function
In this subsection we show that L f (n) has the decreasing increment property with respect to n. That is, L f (n) is the restriction of a concave function to integer values. A by-product of this property is that L f (n) must have a global maximum over the possible values of n. We will call this property "concave" or "concavity" for convenience henceforth. Suppose for now that the possible values of n is a continuous range. First, since R f [j] ∈ {0, · · · , r max }, we have c i ln p Zn (i), (4.11) for some non-negative integer constants c 0 , c 1 , · · · , c rmax such that rmax i=1 c i = k (i.e. the constants represent the empirical pmf of (R f [0], · · · , R f [k − 1])). Therefore a sufficient condition for L f (n) to be concave in n is that ln p Zn (i) is concave in n for each fixed value i ∈ {0, · · · , r max }. Recall the pmf of Z n from (4.8) and (4.9). Let us denote a i F Z (i) and When i = 0, ln p n (i) = ln a i b n i = ln a i + n ln b i , which is linear in n and thus concave in n. When 1 ≤ i ≤ r max , we prove that ln p Zn (i) is concave in n by showing that its second derivative with respect to n is non-positive: We conclude that L f (n) is concave in n.
In Figure 4.2, we show the plot of L f (n) as a function of n for four selected flows whose cardinalities are 150, 611, 1000 and 31536, respectively. The plots are generated based on the same trace data used for We can also observe that in each plot, the value of n at which the curve peaks is close to the actual flow cardinality n f (though not exact the same due to randomness and noise).

Bi-section search implementation
We have shown that L f (n) is concave in n. Since n can only take integer values in {1, 2, · · · , n max }, we can do a bi-section search to find the optimal n.
The pseudo-code above should be self-explanatory. In actual implementation, it is found that the algorithm reaches the optimal n at a faster rate if we replace mid 1 = lb+ub 2 with mid 1 = √ lb · ub, which is equivalent to log mid 1 = 1 2 (log lb + log ub), i.e. a bi-section search on log scale. This faster rate is because, as we already mentioned, most of the flows have small cardinalities. If n max = 10 6 , then at the first while loop iteration, with mid 1 = lb+ub 2 , we have mid 1 = 5 × 10 5 ; but with mid 1 = √ lb · ub we have mid 1 = 1000 -the latter search is quicker if the optimal n is actually small.
With bi-section search, the number of searches required for the estimation of one flow is upper bounded by log n max . This is not ideal compared to the vLL θ estimator which only needs one search/estimation. However, since the estimation process is off-line, the worst-case time complexity of log n max for searching is still acceptable.

Summary of the maximum-likelihood estimator
We wrap up the ideas in the previous sections of this chapter and name this estimator vLL-MLE, meaning that it uses maximum likelihood estimation based on virtual LogLog sketching.
The estimator is summarized as follows. To estimate the cardinality of flow f : 1. Find the values of the k registers in R f : R f [0], · · · , R f [k − 1].
2. From the register values in R, find the CDF of Z by (4.6).
3. Do a bi-section search on integer value n in the range {1, 2, . . . , n max } to maximize the log likelihood function L f (n), output this value of n as the estimate.

Aside: MLE for single-flow estimation
As an aside, we can also use the vLL-MLE estimator for single-flow estimation, with some minor changes. Recall that in single-flow estimation, the goal is to estimate the flow's cardinality n from its k register values. This is similar to the per-flow estimation case, where the goal is to estimate a given flow f 's cardinality n f from its k register values R f [0], · · · , R f [k − 1], except that in the single-flow case, there is no background noise in the register values. That is, the vLL-MLE estimator for per-flow estimation can be used for single-flow estimation by simply letting Z = 0. In this case, we do not need the assumptions related to Z anymore. But the Poissonization approximation of n (i.e. assuming n is large) is still necessary.
To evaluate the performance of the MLE on single-flow estimation, we run experiments to obtain empirical values of its relative standard error for selected values of k and compare it with the HLL estimator. The experiments are performed with n = 10 7 . The results are summarized in Table 4.1. The experiment results show that the MLE's performance is very close to that of HLL, but not any better. It reinforces the claim that the HLL estimator is near-optimal for single-flow estimation.

Experimental performance evaluation
In this section we evaluate the performance of the vLL-MLE estimator and compare it with that of the vLL θ estimator introduced in Chapter 3. Since we have shown that vLL −1 (i.e. vHLL) is the best within the family of the vLL θ estimators, we compare vLL-MLE with vLL −1 in particular. Results for each estimator are generated by running experiments on the same 100 simulated trace files. Figure 4.3 shows the estimation results of both estimators by directly plotting the estimated cardinalities vs. the corresponding actual cardinalities. Each point in the graphs represents one flow, with its x-coordinate value being the actual cardinality of the flow and its y-coordinate value being the estimated cardinality. The more clustered the points are to the equality line y = x, the more accurate the estimator is. We can see that the two estimators have comparable performances.  ) vs. the actual flow cardinality for both estimators. Since there are not many flows for some cardinalities (especial the large cardinalities), we divide the horizontal axis into bins of width 1000 for cardinalities ≤ 10000 and of width 5000 for cardinalities > 10000. In each bin, we calculate the empirical relative bias and relative standard error of the data and interpolate the values for each bin on the graph to form the plots. Again we see that the two estimators have very similar performances.  In both figures, the curves are plotted for cardinalities larger than 1000 only. This is for better graph layouts: we found from the experiment results that estimates for flows with cardinalities less than 1000 are very inaccurate. Accurate estimation for small flows is difficult because of the noise caused by large flows in their register values. Figure 4.4 shows that both estimators are approximately unbiased for large flows (with cardinalities > 1000). Figure 4.5 shows that in general, for both estimators, the estimation accuracy improves as the cardinality gets larger. We do observe a slight curving up at the high end of the graph though, meaning that the accuracy stops improving once the actual flow cardinality reaches a certain level (probably > 40000). Finally, in Figure 4.6 we compare the weighted square error of the vLL-MLE estimator with that of the vLL θ estimators. Experiment results show that the vLL-MLE estimator outperforms the vLL θ estimator for all values of θ; compared to the vLL −1 estimator, the vLL-MLE estimator has a slight improvement of about 3.5%. We conclude that the two estimators, vLL-MLE and vLL −1 (i.e. vHLL), have comparable performances.
identify large flows and estimate their cardinalities directly and more efficiently are useful in this context.