Padua Transforms

Introduction

The Padua transform yields coefficients $a_{ij}$ used to approximate a bivariate function $f(x, y)$ with

\[f(x, y) ≈ \sum_{i, j} a_{ij} \; T_i(y) \; T_j(x)\]

where $T_n(x) := \cos(n \arccos(x))$ is the nth Chebyshev polynomial. For reference, the first few Chebyshev polynomials are

\[\begin{aligned} T_0(x) &= 1 \\ T_1(x) &= x \\ T_2(x) &= 2x^2 - 1 \\ T_3(x) &= 4x^3 - 3x \\ T_4(x) &= 8x^4 - 8x^2 + 1 \end{aligned}\]

The inverse transform takes the coefficients $a_{ij}$ and returns the function $f(x, y)$ evaluated at the so called Padua points.

Start by evaluating a function on the Padua points. Here we approximate the function with a polynomial of total degree 3. The total degree is the degree of the largest polynomial in x plus the degree of the largest polynomial in y.

julia> vals = getpaduapoints(3) do x, y
           y * (2x^2 - 1) + 5 * x * y + 2.5
       end
10-element Vector{Float64}:
  8.5
  2.5
 -3.5
  3.914213562373095
  1.0857864376269049
 -0.5
  2.5
  5.5
 -0.3284271247461903
  5.32842712474619

Then create a PaduaTransformPlan and apply the paduatransform! to get the Chebyshev coefficeints $a_{ij}$.

julia> plan = PaduaTransformPlan{Float64}(3);

julia> coeffs = paduatransform!(zeros(4, 4), plan, vals)
4×4 Matrix{Float64}:
 2.5  0.0  0.0  0.0
 0.0  5.0  1.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

We can go back to values using an InvPaduaTransformPlan and applying invpaduatransform!

julia> invplan = InvPaduaTransformPlan{Float64}(3);

julia> out = invpaduatransform!(Vector{Float64}(undef, getpaduanum(3)), invplan, coeffs)
10-element Vector{Float64}:
  8.5
  2.5
 -3.5
  3.914213562373095
  1.0857864376269049
 -0.5
  2.5
  5.5
 -0.3284271247461903
  5.32842712474619

julia> out ≈ vals
true

where we used getpaduanum to get the number of Padua points and coefficients corresponding to total degree 3.

The Algorithm – Step by Step

To obtain the coefficients $a_{ij}$ we need to evaluate the function $f$ at some points $(x_l, y_k)$ and evaluate

\[a_{ij} = \sum_{k, l} \; f(x_l, y_k) \; T_i(y_k) \; T_j(x_l)\]

If we let $x_l = \cos{\frac{lπ}{n}}$ and $y_k = \cos{\frac{kπ}{m}}$, we have

\[a_{ij} = \sum_{k, l} \; f(\cos{\frac{lπ}{n}}, \cos{\frac{kπ}{m}}) \; T_i(\cos{\frac{kπ}{m}}) \; T_j(\cos{\frac{lπ}{n}})\]

which simplifies to

\[a_{ij} = \sum_{k, l} \; f(\cos{\frac{lπ}{n}}, \cos{\frac{kπ}{m}}) \; \cos{\frac{ikπ}{m}} \; \cos{\frac{jlπ}{n}}\]

because of the definition of the Chebyshev polynomials as $T_n(x) := \cos(n \arccos(x))$. The expression above looks like a discrete cosine transform, which is what we will make use of to implement the Padua transform. Note that this derivation doesn't tell the full story. In fact, we are missing some weighting factors. For further reading on the Padua Transform, please see the following:

For details on the Padua points as good nodes for polynomial interpolation: Marco Caliari, Stefano De Marchi, Marco Vianello. Bivariate polynomial interpolation on the square at new nodal sets

and for details on the implementation of the transform via a discrete cosine transform: Marco Caliari, Stefano De Marchi, Alvise Sommariva, Marco Vianello. Padua2DM: fast interpolation and cubature at the Padua points in Matlab/Octave

The Padua Points

There are multiple ways to define the set of Padua points. For our purposes we will use

\[\textrm{Pad}_n = \{(\cos{\frac{jπ}{n}}, \cos{\frac{iπ}{n + 1}}) \; | \; 0 ≤ j ≤ n, \; 0 ≤ i ≤ n + 1, \; i - j \; \textrm{even} \}\]

where $n$ is the total degree of the polynomial we can approximate using the Padua points. To generate the Padua points we can use the function getpaduapoints as follows

julia> getpaduapoints(3)
10×2 Matrix{Float64}:
  1.0   1.0
  1.0   0.0
  1.0  -1.0
  0.5   0.707107
  0.5  -0.707107
 -0.5   1.0
 -0.5   0.0
 -0.5  -1.0
 -1.0   0.707107
 -1.0  -0.707107

Using do block syntax we can evaluate a function on the Padua points.

julia> vals = getpaduapoints(3) do x, y
           y * (2x^2 - 1) + 5 * x * y + 2.5
       end
10-element Vector{Float64}:
  8.5
  2.5
 -3.5
  3.914213562373095
  1.0857864376269049
 -0.5
  2.5
  5.5
 -0.3284271247461903
  5.32842712474619

The Padua Transform

Having evaluated our function on the Padua points, we can perform the transform. First, we initialise a transform plan.

julia> plan = PaduaTransformPlan{Float64}(3);

The transform consists of writing vals into plan.vals, applying a fast fourier transform and weighting the resulting coefficients to obtain the Chebyshev coefficients. We start by writting values into the values matrix using tovalsmat!.

julia> ChebyshevTransforms.tovalsmat!(plan.vals, vals, 3)
5×4 Matrix{Float64}:
  8.5  0.0      -0.5   0.0
  0.0  3.91421   0.0  -0.328427
  2.5  0.0       2.5   0.0
  0.0  1.08579   0.0   5.32843
 -3.5  0.0       5.5   0.0

vals are written into the matrix plan.vals such that the entry plan.vals[i+1, j+1] corresponds to the Padua point $(\cos{\frac{jπ}{n}}, \cos{\frac{iπ}{n + 1}})$. However, since the Padua points are only those points with $i-j$ even, all entries corresponding to $i-j$ odd are left out. These off grid entries must be filled with 0. Next, we can apply the discrete cosine transform and weight! the coefficients to obtain the Chebyshev coefficients.

julia> plan.dctplan * plan.vals
5×4 Matrix{Float64}:
 60.0   0.0   0.0   0.0
  0.0  30.0   6.0   0.0
  0.0   0.0   0.0   0.0
  0.0   6.0  30.0   0.0
  0.0   0.0   0.0  60.0

julia> ChebyshevTransforms.weight!(plan.vals, 3)
5×4 Matrix{Float64}:
 2.5  0.0  0.0  0.0
 0.0  5.0  1.0  0.0
 0.0  0.0  0.0  0.0
 0.0  1.0  5.0  0.0
 0.0  0.0  0.0  2.5

The weighting factor we apply to the coefficients is

\[w = \frac{1}{n(n+1)} ⋅ \begin{cases} \frac{1}{2} & \textrm{if on vertex} \\ 1 & \textrm{if on edge} \\ 2 & \textrm{if in interior} \\ \end{cases}\]

Finally, we write those coefficients corresponding to a total degree of 3 or lower (the upper left triangular) into the output.

julia> coeffs = zeros(4, 4)
4×4 Matrix{Float64}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

julia> ChebyshevTransforms.fromcoeffsmat!(coeffs, plan.vals, 3)
4×4 Matrix{Float64}:
 2.5  0.0  0.0  0.0
 0.0  5.0  1.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

The Inverse Padua Transform

If we have a set of Chebyshev coefficients and want to obtain the values of the corresponding Chebyshev polynomial on the Padua points, we must use the inverse Padua transform. Again, we start with a plan and write our coefficients into plan.coeffs.

julia> invplan = InvPaduaTransformPlan{Float64}(3);

julia> ChebyshevTransforms.tocoeffsmat!(invplan.coeffs, coeffs)
5×4 Matrix{Float64}:
 2.5  0.0  0.0  0.0
 0.0  5.0  1.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

Next, we weight the coefficients using invweight! and apply a discrete cosine transform.

julia> ChebyshevTransforms.invweight!(invplan.coeffs)
5×4 Matrix{Float64}:
 2.5  0.0   0.0   0.0
 0.0  1.25  0.25  0.0
 0.0  0.0   0.0   0.0
 0.0  0.0   0.0   0.0
 0.0  0.0   0.0   0.0

julia> invplan.dctplan * invplan.coeffs
5×4 Matrix{Float64}:
  8.5      4.5      -0.5      -1.5
  6.74264  3.91421   0.37868  -0.328427
  2.5      2.5       2.5       2.5
 -1.74264  1.08579   4.62132   5.32843
 -3.5      0.5       5.5       6.5

invweight! applies the weighting

\[w = \begin{cases} 1 & \textrm{if on vertex} \\ \frac{1}{2} & \textrm{if on edge} \\ \frac{1}{4} & \textrm{if in interior} \\ \end{cases}\]

Finally, we copy over those values corresponding to $i-j$ even and we have our values back.

julia> out = ChebyshevTransforms.fromvalsmat!(Vector{Float64}(undef, getpaduanum(3)), invplan.coeffs, 3)
10-element Vector{Float64}:
  8.5
  2.5
 -3.5
  3.914213562373095
  1.0857864376269049
 -0.5
  2.5
  5.5
 -0.3284271247461903
  5.32842712474619

julia> vals ≈ out
true

Reference

Padua Points

ChebyshevTransforms.getpaduapoints — Function

getpaduapoints([f::Function,], [T=Float64,] n)

evaluates the function f on the Padua points (of type T) for degree n.

If no function f is provided, getpaduapoints returns a matrix of the paduapoints, where the first and second column represent the x and y components respectively. If f returns a single value, getpaduapoints returns a Vector{T}. If f returns a tuple or other indexable iterable, getpaduapoints returns a Matrix{T} where the i-th column represents the i-th component function of f applied to all of the Padua points.

Examples

julia> getpaduapoints(2)
6×2 Matrix{Float64}:
  1.0   1.0
  1.0  -0.5
  0.0   0.5
  0.0  -1.0
 -1.0   1.0
 -1.0  -0.5

julia> getpaduapoints(Float32, 1)
3×2 Matrix{Float32}:
  1.0   1.0
  1.0  -1.0
 -1.0   0.0

julia> getpaduapoints(Float32, 2) do x, y; x*y; end
6-element Vector{Float32}:
  1.0
 -0.50000006
  0.0
 -0.0
 -1.0
  0.50000006

julia> getpaduapoints(2) do x, y; x*y, 5*x*y; end
6×2 Matrix{Float64}:
  1.0   5.0
 -0.5  -2.5
  0.0   0.0
 -0.0  -0.0
 -1.0  -5.0
  0.5   2.5