Switching Representations via Two-Scale Relations#

The scaling functions at scale \(n\) form an orthonormal basis in function space \(V_n\); these spaces form a hiearchy of nested spaces that converges to the full Hilbert space (\(L^2\)):

\[ V_0 \subset V_1 \subset V_2 \dots \to L^2 \]

The nesting relationship \(V_n \subset V_{n+1}\) is due to the fact that any scaling function \(\phi^n_{il}\) can be written exactly as a linear combination of the scaling functions at scale \(n+1\) whose that are nonzero in its domain, namely \(\phi^{n+1}_{i,2l}\) and \(\phi^{n+1}_{i,2l+1}\). Furthermore, the coefficients in this expansion are scale-indendent, thus the easiest way to illustrate this 2-scale relationship is by expanding the mother scaling functions with the functions at scale 1:

\[ \phi_i(x) = \sqrt{2} \sum_{j=0}^{k-1} \left( h_{ij}^{(0)} \phi^1_{j,0}(x) + h_{ij}^{(1)} \phi^1_{j,1}(x) \right) \]

This should make sense intuitively: functions at scale 1 can represent functions with up to \(k-1\) nodes in each half-interval; the mother polynomials have only up to \((k-1)/2\) nodes in each half-interval, thus they can be represented exactly.

The 2-scale relation between scaling functions can be used to “project” a function represented in terms of scaling functions at scale \(n\) to any scale \(<n\). Clearly any basis transformation from \(V_n\) to \(V_{<n}\) is linear but not invertible; information is lost when projecting onto coarser scales.

It turns out that the information lost when you transform the basis from \(V_1\) to \(V_0\) lives in the function space \(W_0\) spanned by the father wavelets. In other words, \(V_0 \oplus W_0 = V_1\). Similarly, \(V_n \oplus W_n = V_{n+1}\), etc. Since (multi)wavelets are related to the father wavelets by dilations and translations, the 2-scale relationship holds for them also. However, a more important relationship for us is that between scaling functions at finer scale with the wavelets at coarser scale:

\[ \psi_i(x) = \sqrt{2} \sum_{j=0}^{k-1} \left( g_{ij}^{(0)} \phi^1_{j,0}(x) + g_{ij}^{(1)} \phi^1_{j,1}(x) \right) \]

The combination of the two types of two-scale relationships allows to transforms functions represented in \(V_1\) into representation in terms of the bases of \(V_0\) and \(W_0\) exactly, using coefficients that are scale independent. It’s conventional to express such transformation in a matrix form:

\[\begin{split} \pmatrix{\boldsymbol{\phi}^n_{2l}(x) \\ \boldsymbol{\phi}^n_{2l+1}(x)} = \pmatrix{\boldsymbol{H}^{(0)} & \boldsymbol{G}^{(0)} \\ \boldsymbol{H}^{(1)} & \boldsymbol{G}^{(1)}} \pmatrix{\boldsymbol{\phi}^{n-1}_{l}(x) \\ \boldsymbol{\psi}^{n-1}_{l}(x)} \end{split}\]

where we have used bold symbols to the column vector of \(k\) scaling or wavelet functions in a given node.

In the code cell below you will find a function which will return a selected filter (\(\boldsymbol{H}^{(i)}\) or \(\boldsymbol{G}^{(i)}\)) for a given basis (Legendre/Interpolating), a given order (limited to \(k \le 10\)).

import sys, struct
import numpy as np

# Gets one of the four filter matrices
# Basis: L (legendre) or I (interpolating)
# Type: H or G
# Flag: 0 or 1
# n: from 1 to 11
def getfilter(basis = "L", ftype = "H", flag = 0, n = 4):
    order = n - 1
    name = "mwfilters/" + basis + "_" + ftype + "0_" + str(order)
    filterFile = open(name, 'rb')
    filterData = filterFile.read()
    filterList = []
    for i in range((n)*(n)):
        rstart = i * 8
        rend = (i + 1) * 8
        filterList.append(struct.unpack_from("d",filterData[rstart:rend]))
    filterArray = np.array(filterList)
    filterMatrix = np.ndarray((n, n),buffer=filterArray)
    if (flag == 1):
        if (ftype == "G" and n%2 != 0 ):
            sign_np = -1
        else:
            sign_np = 1
        for i in range(n):
            for j in range(n):
                sign_ij = 1
                if (((i+j)%2) == 0):
                    sign_ij = -1
                sign = sign_ij * sign_np
                filterMatrix[i,j] *= sign
    return filterMatrix

#filters for a third order (k+1=4) Legendre basis

basis = "L" # alternative "I" for interpolating
n = 4 # n = k+1  (avaliable up to k=10)
H0 = getfilter(basis, "H", 0, n)
H1 = getfilter(basis, "H", 1, n)
G0 = getfilter(basis, "G", 0, n)
G1 = getfilter(basis, "G", 1, n)
print(H0)
print(H1)
print(G0)
print(G1)

[[ 7.07106781e-01  0.00000000e+00  4.88535468e-77  0.00000000e+00]
 [-6.12372436e-01  3.53553391e-01 -4.48460293e-77  2.23275976e-77]
 [ 9.16004002e-78 -6.84653197e-01  1.76776695e-01 -3.96935068e-77]
 [ 2.33853587e-01  4.05046294e-01 -5.22912517e-01  8.83883476e-02]]
[[-7.07106781e-01  0.00000000e+00 -4.88535468e-77  0.00000000e+00]
 [-6.12372436e-01 -3.53553391e-01 -4.48460293e-77 -2.23275976e-77]
 [-9.16004002e-78 -6.84653197e-01 -1.76776695e-01 -3.96935068e-77]
 [ 2.33853587e-01 -4.05046294e-01 -5.22912517e-01 -8.83883476e-02]]
[[ 6.54942862e-76 -1.53392998e-01 -5.94088526e-01  3.51467512e-01]
 [-1.54303350e-01 -2.67261242e-01 -1.72516390e-01  6.12372436e-01]
 [-1.98162199e-75 -8.78668779e-02 -3.40306955e-01 -6.13571991e-01]
 [ 2.15645487e-01  3.73508940e-01  4.43622131e-01  3.42326598e-01]]
[[-6.54942862e-76 -1.53392998e-01  5.94088526e-01  3.51467512e-01]
 [-1.54303350e-01  2.67261242e-01 -1.72516390e-01 -6.12372436e-01]
 [ 1.98162199e-75 -8.78668779e-02  3.40306955e-01 -6.13571991e-01]
 [ 2.15645487e-01 -3.73508940e-01  4.43622131e-01 -3.42326598e-01]]

Exercises#

  • Show (numerically) that the four filter matrices constitute a unitary transformation.

  • Use the dilation and translation relationships to obtain, for a given scaling/wavelet function, the corresponding function at the next finer scale (larger \(n\)).

  • Use the two-scale relationship to verify that the scaling function obtained by translation/dilation can also be obtained by linear transformation of the scaling&wavelet functions at the next coarser scale (smaller \(n\)).

from vampyr import LegendreBasis, InterpolatingBasis
import matplotlib.pyplot as plt
import numpy as np

# Create Legendre basis of order 4 and Interpolating basis of order 3
legendre_4 = LegendreBasis(order=4)
interpol_3 = InterpolatingBasis(order=3)