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.