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

Numbers of Coefficients and Points

ChebyshevTransforms.getpaduanum — Function

getpaduanum(n)

calculates number of Padua points needed to approximate a function using Chebyshev polynomials up to total degree n. This number is equal to the number of coefficients. The formula is

\[N = (n + 1) ⋅ (n + 2) ÷ 2\]

Examples

julia> getpaduanum(13)
105

ChebyshevTransforms.getdegree — Function

getdegree(N)

calculates total degree, given the number of coefficients or Padua points N. Throws an error if N is not a possible number of Padua points. The formula is

\[n = \frac{\sqrt{1 + 8N} - 3}{2}\]

Examples

julia> getdegree(105)
13

julia> getdegree(104)
ERROR: ArgumentError: number of Padua points or coeffs must be (n + 1) * (n + 2) ÷ 2
[...]

ChebyshevTransforms.nextpaduanum — Function

nextpaduanum(N)

get next valid number of Padua points ≥ N.

Examples

julia> nextpaduanum(104)
105

ChebyshevTransforms.nextdegree — Function

nextdegree(N)

get degree of nextpaduanum(N).

Examples

julia> nextdegree(104)
13

Padua Transform

ChebyshevTransforms.PaduaTransformPlan — Type

PaduaTransformPlan{T}(n::Integer)

create plan to compute coefficients of Chebyshev polynomials in 2D up to total degree n using the Padua transform.

ChebyshevTransforms.paduatransform! — Function

paduatransform!(coeffs, P::PaduaTransformPlan, vals[, lex])

obtain coefficients of Chebyshev polynomials on 2D via the Padua transform, given values vals evaluated at the Padua points. Coefficients will be written into coeffs, which should either be a matrix, a vector or an iterable of either.

If an iterable of vals vectors and a corresponding iterable of coeffs matrixes or vectors is given, each vals vector will be transformed seperately in a multidimensional Padua transform.

lex determines the order in which coefficients are written into out if out is a vector. See fromcoeffsmat! for details.

Warning

if coeffs is a matrix, make sure that all entries in the lower right diagonal are zero as these will not get overwritten.

Examples

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

julia> f(x, y) = 3 + 4x + 5 * x*y, 6 + 7y
f (generic function with 1 method)

julia> vals = getpaduapoints(f, 2)
6×2 Matrix{Float64}:
 12.0  13.0
  4.5   2.5
  3.0   9.5
  3.0  -1.0
 -6.0  13.0
  1.5   2.5

julia> paduatransform!(zeros(3, 3), plan, vals[:, 1])
3×3 Matrix{Float64}:
 3.0  4.0  0.0
 0.0  5.0  0.0
 0.0  0.0  0.0

julia> paduatransform!(zeros(6), plan, vals[:, 2], Val(true))
6-element Vector{Float64}:
 6.0
 0.0
 7.0
 0.0
 0.0
 0.0

julia> paduatransform!((zeros(3, 3), zeros(3, 3)), plan, vals)
([3.0 4.0 0.0; 0.0 5.0 0.0; 0.0 0.0 0.0], [6.0 0.0 0.0; 7.0 0.0 0.0; 0.0 0.0 0.0])

Inverse Padua Transform

ChebyshevTransforms.InvPaduaTransformPlan — Type

InvPaduaTransformPlan{T}(n::Integer)

create plan to compute values on Padua points, given coefficients of Chebyshev polynomials up to total degree n.

ChebyshevTransforms.invpaduatransform! — Function

invpaduatransform!(vals::AbstractVector, IP::InvPaduaTransformPlan, coeffs::AbstractMatrix)

evaluates the polynomial defined by the coefficients of Chebyshev polynomials coeffs on the Padua points using the inverse transform plan IP and writes the resulting values into vals.

If an iterable of vals vectors or a vals matrix and a corresponding iterable of coeffs matrixes is given, each coeffs matrix will be transformed seperately.

Examples

julia> iplan = InvPaduaTransformPlan{Float64}(2);

julia> coeffs = [[3 4 0; 0 5 0; 0 0 0], [3 0 1; 4 0 0; 0 0 0]]
2-element Vector{Matrix{Int64}}:
 [3 4 0; 0 5 0; 0 0 0]
 [3 0 1; 4 0 0; 0 0 0]

julia> invpaduatransform!(zeros(6), iplan, coeffs[1])
6-element Vector{Float64}:
 12.0
  4.5
  3.0
  3.0
 -6.0
  1.5

julia> invpaduatransform!((zeros(6), zeros(6)), iplan, coeffs)
([12.0, 4.5, 3.0, 3.0, -6.0, 1.5], [8.0, 2.0, 4.0, -2.0, 8.0, 2.0])

julia> invpaduatransform!(zeros(6, 2), iplan, coeffs)
6×2 Matrix{Float64}:
 12.0   8.0
  4.5   2.0
  3.0   4.0
  3.0  -2.0
 -6.0   8.0
  1.5   2.0

julia> getpaduapoints(2) do x, y; (3 + 4x + 5 * x*y, 3 + 2x^2 - 1 + 4y); end
6×2 Matrix{Float64}:
 12.0   8.0
  4.5   2.0
  3.0   4.0
  3.0  -2.0
 -6.0   8.0
  1.5   2.0

Internals

ChebyshevTransforms.paduapoint — Function

paduapoint(T::Type, j::Integer, i::Integer, n::Integer)

returns the Padua point $z_{ij}$, where

\[z_{ij} = (\cos{\frac{jπ}{n}}, \cos{\frac{iπ}{n+1}})\]

Note, that only points with $i-j$ even are actually Padua points. Check with ispadua.

Examples

julia> [ChebyshevTransforms.paduapoint(Float32, x, y, 1) for y in 0:1+1, x in 0:1]
3×2 Matrix{Tuple{Float32, Float32}}:
 (1.0, 1.0)   (-1.0, 1.0)
 (1.0, 0.0)   (-1.0, 0.0)
 (1.0, -1.0)  (-1.0, -1.0)

ChebyshevTransforms.ispadua — Function

ispadua(i, j)

returns if paduapoint at position (i, j) is a Padua point.

Examples

julia> pointornothing(i, j, n) = ChebyshevTransforms.ispadua(i, j) ? ChebyshevTransforms.paduapoint(Float64, j, i, n) : nothing
pointornothing (generic function with 1 method)

julia> [pointornothing(y, x, 2) for y in 0:3, x in 0:2]
4×3 Matrix{Union{Nothing, Tuple{Float64, Float64}}}:
 (1.0, 1.0)   nothing      (-1.0, 1.0)
 nothing      (0.0, 0.5)   nothing
 (1.0, -0.5)  nothing      (-1.0, -0.5)
 nothing      (0.0, -1.0)  nothing

ChebyshevTransforms.tovalsmat! — Function

tovalsmat!(mat::Matrix, from::AbstractVector, degree::Integer)

write values of function evaluated at Padua points from from to matrix mat.

Examples

julia> ChebyshevTransforms.tovalsmat!(ones(3 + 2, 3 + 1), 1:getpaduanum(3), 3)
5×4 Matrix{Float64}:
 1.0  0.0  6.0   0.0
 0.0  4.0  0.0   9.0
 2.0  0.0  7.0   0.0
 0.0  5.0  0.0  10.0
 3.0  0.0  8.0   0.0

julia> ChebyshevTransforms.tovalsmat!(ones(2 + 2, 2 + 1), 1:getpaduanum(2), 2)
4×3 Matrix{Float64}:
 1.0  0.0  5.0
 0.0  3.0  0.0
 2.0  0.0  6.0
 0.0  4.0  0.0

ChebyshevTransforms.weight! — Function

weight!(mat::AbstractMatrix, degree::Integer)

weight fourier coefficients to obtain Chebyshev coefficients as part of a paduatransform!. The weighting factor applied 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}\]

Examples

julia> weight!(ones(4+2, 4+1), 4)
6×5 Matrix{Float64}:
 0.025  0.05  0.05  0.05  0.025
 0.05   0.1   0.1   0.1   0.05
 0.05   0.1   0.1   0.1   0.05
 0.05   0.1   0.1   0.1   0.05
 0.05   0.1   0.1   0.1   0.05
 0.025  0.05  0.05  0.05  0.025

ChebyshevTransforms.fromcoeffsmat! — Function

fromcoeffsmat!(to::AbstractVector, mat::Matrix, degree::Integer, ::Val{lex})

write Chebyshev coefficients from mat into vector to. lex::Bool determines whether coefficients should be written in lexigraphical order or not. The lower right triangle does not get written into to. These would represent greater polynomial degrees than degree.

If lex is Val(true) the coefficients correspond to the following basis polynomials

\[T_0(x) T_0(y), T_1(x) T_0(y), T_0(x) T_1(y), T_2(x) T_0(y), T_1(x) T_1(y), T_0(x) T_2(y), ...\]

else if lex is Val(false) they correspond to

\[T_0(x) T_0(y), T_0(x) T_1(y), T_1(x) T_0(y), T_0(x) T_2(y), T_1(x) T_1(y), T_2(x) T_0(y), ...\]

Examples

julia> mat = [(x, y) for y in 0:2+1, x in 0:2]
4×3 Matrix{Tuple{Int64, Int64}}:
 (0, 0)  (1, 0)  (2, 0)
 (0, 1)  (1, 1)  (2, 1)
 (0, 2)  (1, 2)  (2, 2)
 (0, 3)  (1, 3)  (2, 3)

julia> to1 = similar(mat, getpaduanum(2)); to2 = similar(mat, getpaduanum(2));

julia> ChebyshevTransforms.fromcoeffsmat!(to1, mat, 2, Val(true))
6-element Vector{Tuple{Int64, Int64}}:
 (0, 0)
 (1, 0)
 (0, 1)
 (2, 0)
 (1, 1)
 (0, 2)

julia> ChebyshevTransforms.fromcoeffsmat!(to2, mat, 2, Val(false))
6-element Vector{Tuple{Int64, Int64}}:
 (0, 0)
 (0, 1)
 (1, 0)
 (0, 2)
 (1, 1)
 (2, 0)

fromcoeffsmat!(to::AbstractMatrix, mat::Matrix, degree::Integer)

copy Chebyshev coefficients from mat to to without copying coefficients coresponding to total degree greater than degree.

Examples

julia> ChebyshevTransforms.fromcoeffsmat!(zeros(4, 4), reshape(1:20, 5, 4), 3)
4×4 Matrix{Float64}:
 1.0  6.0  11.0  16.0
 2.0  7.0  12.0   0.0
 3.0  8.0   0.0   0.0
 4.0  0.0   0.0   0.0

ChebyshevTransforms.tocoeffsmat! — Function

tocoeffsmat!(mat::AbstractMatrix, coeffs::AbstractMatrix)

writes coefficients in coeffs into matrix mat for the invpaduatransform!.

Examples

julia> ChebyshevTransforms.tocoeffsmat!(zeros(5, 4), reshape(1:16, 4, 4))
5×4 Matrix{Float64}:
 1.0  5.0   9.0  13.0
 2.0  6.0  10.0  14.0
 3.0  7.0  11.0  15.0
 4.0  8.0  12.0  16.0
 0.0  0.0   0.0   0.0

ChebyshevTransforms.invweight! — Function

invweight!(coeffs::AbstractMatrix)

weight Chebyshev coefficients before the Fourier transform for the invpaduatransform!. using 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}\]

Examples

julia> ChebyshevTransforms.invweight!(ones(5, 5))
5×5 Matrix{Float64}:
 1.0  0.5   0.5   0.5   1.0
 0.5  0.25  0.25  0.25  0.5
 0.5  0.25  0.25  0.25  0.5
 0.5  0.25  0.25  0.25  0.5
 1.0  0.5   0.5   0.5   1.0

ChebyshevTransforms.fromvalsmat! — Function

fromvalsmat!(to::AbstractVector, mat::AbstractMatrix, n::Integer)

write values from mat into the vector to after an invpaduatransform! of total degree n.

Examples

julia> ChebyshevTransforms.fromvalsmat!(zeros(10), reshape(1:20, 5, 4), 3)
10-element Vector{Float64}:
  1.0
  3.0
  5.0
  7.0
  9.0
 11.0
 13.0
 15.0
 17.0
 19.0