<!-- dom:TITLE: Demo - Some fast transforms -->
# Demo - Some fast transforms
<!-- dom:AUTHOR: Mikael Mortensen Email:mikaem@math.uio.no at Department of Mathematics, University of Oslo. -->
<!-- Author: -->  
**Mikael Mortensen** (email: `mikaem@math.uio.no`), Department of Mathematics, University of Oslo.

Date: **May 27, 2021**

**Summary.** This demo will show how to compute fast forward transforms for the three
different Dirichlet bases that are implemented for Chebyshev
polynomials in Shenfun.












## Forward and backward transforms

A function $u(x)$ can be approximated in a finite global spectral
expansion $u_N(x)$ as

<!-- Equation labels as ordinary links -->
<a id="eq:expansion"></a>

$$
\label{eq:expansion} \tag{1}
u_N(x) = \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x), \quad \forall \, x \, \in [-1, 1],
$$

where $\phi_k(x)$ are the basis functions and $\boldsymbol{\hat{u}} = \{\hat{u}_k\}_{k=0}^{N-1}$
are the expansion coefficients. The function $u_N(x)$ is continuous
on the interval domain $[-1, 1]$. The span of the basis functions
$V_N = \text{span} \{\phi_k\}_{k=0}^{N-1}$ represents a functionspace.
Associated with this functionspace is a set of quadrature points
$\{x_k\}_{k=0}^{N-1}$ that, along with quadrature weights $\{\omega_k\}_{k=0}^{N-1}$, can be used
for efficient integration. We can also evaluate the function $u_N(x)$ at
these quadrature points to get the sequence
$\boldsymbol{u} = \{u_N(x_k)\}_{k=0}^{N-1}$. If $\boldsymbol{\hat{u}}=\{\hat{u}_k\}_{k=0}^{N-1}$ are *known*,
then $\boldsymbol{u}$ can be evaluated directly from
Eq. ([1](#eq:expansion))

<!-- Equation labels as ordinary links -->
<a id="eq:expansionQ"></a>

$$
\label{eq:expansionQ} \tag{2}
u_N(x_j) = \sum_{k=0}^{N-1} \hat{u}_k \phi_k(x_j), \quad \forall \, j=0,1, \ldots, N-1.
$$

This would correspond to a *backward* transform according to
the Shenfun terminology. A direct evaluation of the backward
([2](#eq:expansionQ)) transform takes $\mathcal{O}(N^2)$
operations since it requires a double sum (over both $j$
and $k$). A *fast* transform is
a transform that can be computed in $\mathcal{O}(N \log N)$ operations.
This is what the Fast Fourier Transform (FFT) does. It computes a double
sum, like ([2](#eq:expansionQ)), in $\mathcal{O}(N \log N)$ operations.


The other way around, computing $\{\hat{u}_k\}_{k=0}^{N-1}$ from the
*known* $\{u_N(x_k)\}_{k=0}^{N-1}$ corresponds to a *forward* transform.
The forward transform is computed using a projection of $u$
into $V_N$, which is formulated as: find $u_N \in V_N$ such that

<!-- Equation labels as ordinary links -->
<a id="eq:projection"></a>

$$
\label{eq:projection} \tag{3}
(u_N-u, v)_{\omega^{\sigma}} = 0, \quad \forall \, v \in V_{N},
$$

where $(a, b)_{\omega^{\sigma}} = \int_{I} a b \omega^{\sigma} dx$ is the
inner product in $L^2_{\omega^{\sigma}}(I)$, and $\omega^{\sigma}(x)=(1-x^2)^{\sigma}$  is a weight function.
For Chebyshev polynomials the weight function is usually $\omega^{-1/2}=(1-x^2)^{-1/2}$.
Inserting for $u_N$ and $v=\phi_k$, we get

<!-- Equation labels as ordinary links -->
<a id="_auto1"></a>

$$
\begin{equation}
\sum_{j=0}^{N-1}(\phi_j, \phi_k)_{\omega^{\sigma}} \hat{u}_{j} = (u, \phi_k)_{\omega^{\sigma}}, 
\label{_auto1} \tag{4}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<a id="_auto2"></a>

$$
\begin{equation}  
B \boldsymbol{\hat{u}} = \boldsymbol{\tilde{u}}, 
\label{_auto2} \tag{5}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<a id="_auto3"></a>

$$
\begin{equation}  
\boldsymbol{\hat{u}} = B^{-1} \boldsymbol{\tilde{u}},
\label{_auto3} \tag{6}
\end{equation}
$$

where
$\boldsymbol{\tilde{u}} = \{(u, \phi_k)_{\omega^{\sigma}}\}_{k=0}^{N-1}$ and the mass matrix
$B = (b_{kj})_{k,j=0}^{N-1}$, with $b_{kj}=(\phi_j, \phi_k)_{\omega^{\sigma}}$.

Note that the forward transform requires both an inner product
$\boldsymbol{\tilde{u}}$ and a matrix inversion. By a *fast* forward transform
we mean a transform that can be computed in $\mathcal{O}(N \log N)$
operations. If $B$ is a diagonal or banded matrix, the matrix inversion costs $\mathcal{O}(N)$,
and the limiting factor is then the inner product. Like for the backward transform,
the inner product, computed with quadrature, is a double sum

$$
(u, \phi_k)_{\omega^{\sigma}} = \sum_{j=0}^{N-1} u(x_j) \phi_k(x_j) \omega_j, \quad \forall \, k = 0, 1, \ldots, N-1,
$$

where $\{\omega_j\}_{j=0}^{N-1}$ are the quadrature weights.

A naive implementation of the inner product
takes $\mathcal{O}(N^2)$ operations. However,
for Chebyshev polynomials we can compute the double loop with
fast $\mathcal{O}(N \log N)$ discrete sine or cosine transforms,
that are versions of the FFT. To see this, assume that the basis functions are $\phi_k(x) =T_k(x)$, where
$T_k(x)$ is the $k$'th Chebyshev polynomial of the first kind,
and the weight function is $\omega^{-1/2}$.
We then choose Gauss-Chebyshev points $x_j = \cos(\theta_j)$,
where $\theta_j=\pi (2j+1)/(2N)$, and the associated quadrature weights
that are constant $\omega_j = \pi/N$. The Chebyshev polynomials evaluated
on the quadrature points can now
alternatively be written as $T_k(x_j) = \cos(k \theta_j)$,
and the inner product becomes

$$
(u, T_k)_{\omega^{-1/2}} = \sum_{j=0}^{N-1} u(x_j) \cos(k \theta_j) \pi/N, \quad \forall \, k = 0, 1, \ldots, N-1.
$$

From the [FFTW documentation](http://www.fftw.org/fftw3_doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html#g_t1d-Real_002deven-DFTs-_0028DCTs_0029)
we recognise this sum as half a DCT-II (the FFTW DCT-II has a factor
2 in front of the sum) of $\boldsymbol{u}\pi/N$. Hence, we can compute the inner product as

$$
(u, T_k)_{\omega^{-1/2}} = \frac{\pi}{2N} \text{dct}^{II}(\boldsymbol{u})_k, \quad k = 0, 1, \ldots, N-1.
$$

## Dirichlet bases

The basis function $T_k$ satisfies $T_k(\pm 1) = (\pm 1)^k$ at the
boundaries of the domain, and the space $S_N=\text{span}\{T_k\}_{k=0}^{N-1}$,
of dimension $N$,
is thus not associated with any specific set of boundary conditions.
A functionspace for homogeneous Dirichlet boundary conditions is
given as $V_N=\{v\in S_N | v(\pm 1)=0 \}$. Because of the two restrictions
the space has dimension $N-2$.

There are several different choices of basis functions
for $V_N$.
The most interesting we name $\phi_k^n$, for integer $n$, and
define them as

<!-- Equation labels as ordinary links -->
<a id="_auto4"></a>

$$
\begin{equation}
\phi^n_k = \omega T^{(n)}_{k+n} = (1-x^2) T^{(n)}_{k+n},
\label{_auto4} \tag{7}
\end{equation}
$$

where $T^{(n)}_{k+n}$ is the $n$'th derivative of $T_{k+n}$. We have
for any integer $n$ that $V_N=\text{span}\{\phi^n_k\}_{k=0}^{N-3}$, and an
expansion in any of these basis functions is

<!-- Equation labels as ordinary links -->
<a id="eq:uNgeneric"></a>

$$
\begin{equation}
\label{eq:uNgeneric} \tag{8}
    u_N = \sum_{k=0}^{N-3} \hat{u}^n_k \phi^n_k.
\end{equation}
$$

We can find the sequence $\{\hat{u}^n_{k}\}_{k=0}^{N-3}$ for any $n$
using a projection into the space $V_N$. The projection is computed
by using Eq. ([8](#eq:uNgeneric)) and $v=\phi^n_k$ in
Eq. ([3](#eq:projection))

<!-- Equation labels as ordinary links -->
<a id="eq:projortho"></a>

$$
\begin{equation}
\label{eq:projortho} \tag{9}
   \sum_{j=0}^{N-3} ( T^{(n)}_{j+n},  T^{(n)}_{k+n})_{\omega^{\sigma+2}} \hat{u}^{n}_j = (u, T^{(n)}_{k+n})_{\omega^{\sigma+1}}.
\end{equation}
$$

Now how can this projection be computed as efficiently as possible?
The Chebyshev polynomials and their derivatives are known to satisfy
the following orthogonality relation

<!-- Equation labels as ordinary links -->
<a id="eq:orthon"></a>

$$
\begin{equation}
\label{eq:orthon} \tag{10}
  (T^{(n)}_j, T^{(n)}_k)_{\omega^{n-1/2}} = \alpha^{n}_k \delta_{kj}, \quad \text{for}\, n \ge 0,
\end{equation}
$$

where $\delta_{kj}$ is the Kronecker delta function and

<!-- Equation labels as ordinary links -->
<a id="_auto5"></a>

$$
\begin{equation}
    \alpha^n_k = \frac{c_{k+n}\pi k (k+n-1)!}{2(k-n)!},
\label{_auto5} \tag{11}
\end{equation}
$$

where $c_0=2$ and $c_k=1$ for $k>0$. This can be used in
computing ([9](#eq:projortho)), because we just
need to choose the $\sigma$ that leads to a diagonal mass matrix.
For $n=(0, 1, 2)$ this will be $\sigma=-5/2, -3/2$ and $-1/2$,
respectively. So, choosing $\sigma=-5/2, -3/2$ and $-1/2$
for $n=0, 1$ and 2, respectively, will lead to a diagonal
mass matrix $( T^{(n)}_{j+n},  T^{(n)}_{k+n})_{\omega^{\sigma+2}}$.
Using these $\sigma$'s we can invert the diagonal mass matrices
in Eq. ([9](#eq:projortho)) to get

<!-- Equation labels as ordinary links -->
<a id="_auto6"></a>

$$
\begin{equation}
    \hat{u}^n_k = \frac{1}{\alpha^n_{k+n}}(u, T^{(n)}_{k+n})_{\omega^{\sigma+1}}, \quad k=0, 1, \ldots, N-3, \text{ for } n \in (0, 1, 2).
\label{_auto6} \tag{12}
\end{equation}
$$

Using now quadrature, $1-x^2_i=\sin^2 \theta_i$ and the
fast transforms $(u, T_k)_{\omega^{-1/2}} = \pi/2/N \text{dct}^{II}(\boldsymbol{u})_k$
and  $(u, U_k)_{\omega^{-1/2}} = \pi/2/N \text{dst}^{II}(\boldsymbol{u}/\sin \boldsymbol{\theta})_k$,
where $\boldsymbol{u}/\sin \boldsymbol{\theta}$ implies element-wise division,
we get

<!-- Equation labels as ordinary links -->
<a id="eq:fast1"></a>

$$
\begin{equation}
    \hat{u}^0_k = \frac{1}{c_k N} \text{dct}^{II}(\boldsymbol{u}/\sin^2 \boldsymbol{\theta})_k, \quad k = 0, 1, \ldots, N-3, \label{eq:fast1} \tag{13} 
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<a id="eq:fast2"></a>

$$
\begin{equation}  
    \hat{u}^1_k = \frac{1}{(k+1)N}\text{dst}^{II}(\boldsymbol{u}/\sin \boldsymbol{\theta})_k, \quad k = 0, 1, \ldots, N-3, \label{eq:fast2} \tag{14}
\end{equation}
$$

<!-- Equation labels as ordinary links -->
<a id="eq:fast3"></a>

$$
\begin{equation}  
    \hat{u}^2_k = \frac{1}{2(k+2)}\left(\hat{u}^1_k - \hat{u}^1_{k+2} \right), \quad k=0, 1, \ldots, N-3. \label{eq:fast3} \tag{15}
\end{equation}
$$

The last one requires some work, using the identity
$\phi^2_k=(1-x^2)T''_{k+2}=0.5(k+2)(k+3)(U_k - (k+1)/(k+3)U_{k+2})$.

## Verification

To validate all the fast methods we compute the projection first regularly
using the Shenfun function [project](https://github.com/spectralDNS/shenfun/blob/master/shenfun/forms/project.py),
which is using $\sigma=-1/2$, and then the fast methods above. The two
projections should be the same, but they will not give identical results.
In general, the fast transforms above should be both faster and more
accurate, because they only take a discrete transform and merely a diagonal
mass matrix inversion.

We start the implementation by importing necessary modules from Shenfun
and mpi4py-fft


```
from shenfun import *
from mpi4py_fft import fftw
```

The three bases $\{\phi^n_k\}_{k=0}^{N-3}$ are implemented
with slightly different scaling in shenfun.
The first, with $n=0$, is obtained with no special scaling using


```
N = 20
D0 = FunctionSpace(N, 'C', bc=(0, 0), basis='Heinrichs')
```

The second basis is implemented in Shenfun as $\phi_k = \frac{2}{k+1}\phi^1_k$,
which can be simplified as

<!-- Equation labels as ordinary links -->
<a id="eq:ft:shen"></a>

$$
\label{eq:ft:shen} \tag{16}
\phi_k(x) = T_k-T_{k+2}, \quad k=0,1, \ldots, N-3,
$$

and implemented as


```
D1 = FunctionSpace(N, 'C', bc=(0, 0)) # this is the default basis
```

Because of the scaling the expansion coefficients for $\phi_k$  are
$\hat{u}^{\phi}_k=\frac{k+1}{2}\hat{u}^1_k$. Using ([14](#eq:fast2)) we get

$$
\hat{u}^{\phi}_k = \frac{1}{2N}\text{dst}^{II}(\boldsymbol{u}/\sin \boldsymbol{\theta})_k, \quad k = 0, 1, \ldots, N-3.
$$

The third basis is also scaled and implemented in Shenfun as $\psi_k = \frac{2}{(k+3)(k+2)}\phi^2_k$,
which can be simplified using Chebyshev polynomials of the second
kind $U_k$

<!-- Equation labels as ordinary links -->
<a id="eq:ft:dirichletU"></a>

$$
\label{eq:ft:dirichletU} \tag{17}
\psi_k(x) = U_k-\frac{k+1}{k+3}U_{k+2}, \quad k=0,1, \ldots, N-3.
$$

We get the basis using


```
D2 = FunctionSpace(N, 'C', bc=(0, 0), basis='DirichletU')
```

and the expansion coefficients are found as
$\hat{u}^{\psi}_k = \frac{(k+3)(k+2)}{2} \hat{u}^2_k$.

For verification of all the fast transforms we first create a vector
consisting of random expansion coefficients, and then transform
it backwards to physical space


```
f = Function(D0, buffer=np.random.random(N))
f[-2:] = 0
fb = f.backward().copy()
```

Next, we perform the regular projections into the three spaces
`D0`, `D1` and `D2`, using the inner product
in $L^2_{\omega^{-1/2}}$. Now `u0`, `u1` and `u2` will be the
three solution vectors
$\boldsymbol{\hat{u}}^{\varphi}$, $\boldsymbol{\hat{u}}^{\phi}$
and $\boldsymbol{\hat{u}}^{\psi}$, respectively.


```
u0 = project(fb, D0)
u1 = project(fb, D1)
u2 = project(fb, D2)
```

Now compute the fast transforms and assert that they are equal to `u0`, `u1` and `u2`


```
theta = np.pi*(2*np.arange(N)+1)/(2*N)
# Test for n=0
dct = fftw.dctn(fb.copy(), type=2)
ck = np.ones(N); ck[0] = 2
d0 = dct(fb/np.sin(theta)**2)/(ck*N)
assert np.linalg.norm(d0-u0) < 1e-8, np.linalg.norm(d0-f0)
# Test for n=1
dst = fftw.dstn(fb.copy(), type=2)
d1 = dst(fb/np.sin(theta))/(2*N)
assert np.linalg.norm(d1-u1) < 1e-8
# Test for n=2
ut = d1
k = np.arange(N)
d2 = Function(D2)
d2[:-2] = (k[:-2]+3)/2/(k[:-2]+1)*ut[:-2]
d2[:-2] = d2[:-2] - 0.5*ut[2:]
assert np.linalg.norm(d2-u2) < 1e-8
```

That's it! If you make it to here with no errors, then the three tests pass, and the fast transforms are equal to the slow ones, at least within given precision.

Let's try some timings


```
%timeit project(fb, D1)
```


```
%timeit dst(fb/np.sin(theta))/(2*N)
```

We can precompute the sine term, because it does not change


```
dd = np.sin(theta)*2*N
%timeit dst(fb/dd)
```

The other two transforms are approximately the same speed.


```
%timeit dct(fb/np.sin(theta)**2)/(ck*N)
```

<!-- ======= Bibliography ======= -->

f7ea705bd494e461b85a50b43d83ff6653c93da0
ipynb
Jupyter Notebook
content/fasttransforms.ipynb
mikaem/shenfun-demos
c2ad13d62866e0812068673fdb6a7ef68ecfb7f2
BSD-2-Clause
content/fasttransforms.ipynb
mikaem/shenfun-demos
c2ad13d62866e0812068673fdb6a7ef68ecfb7f2
BSD-2-Clause
2021-09-21T16:10:01.000Z
2021-09-21T16:10:01.000Z
content/fasttransforms.ipynb
mikaem/shenfun-demos
c2ad13d62866e0812068673fdb6a7ef68ecfb7f2
BSD-2-Clause
Qwen/Qwen-72B
1. YES
2. YES