Multiresolution regression of a function with non-continuous derivative

Author: Daniele E. Schiavazzi

Date: June 28th, 2017

References

• Schiavazzi D., Doostan A., Iaccarino G. and Marsden A.L., A Generalized Multi-resolution Expansion for Uncertainty Propagation with Application to Cardiovascular Modeling, Computer Methods in Applied Mechanics and Engineering, 314(1):196-221, 2017. Link
• Schiavazzi D., Doostan A. and Iaccarino G., Sparse multiresolution regression for uncertainty propagation, International Journal for Uncertainty Quantification, 4(4):303-331, 2014. Link

Objectives of this tutorial:

• Introduce multiwavelet basis for stochastic regression.
• Show how to play with the parameters needed to generate these functions.
• Demonstrate application in regression for a simple 1D function with discontinuous derivative.

Let's first define a simple function $$ f(x) = \frac{1.0}{\vert 0.25 - x^2\vert + 0.1} $$ defined in the $[0,1]$ interval, with discontinuous slope at 0.5. Let's plot it to see what it looks like.

Define a flag to switch between random and equally spaced sampling points. This allows to quickly check what happens using random instead of equally spaced points and by changing their number.

Let's visualize the selected sampling points

Since we would like to perform regression with Legendre polynomials first, we need to scale xVals to fit in the $[-1,1]$ interval.

We want to perform regression using Relevance Vector Machine regression. So once again we initialize the associated python class as for the previous tutorial.

Let's increase the order of the polynomials regression to 30.

Rudiments of Multiresolution analysis and multiwavelet basis

Multiresolution analysis

A multiresolution approximation of $\mathbf{L}^2([0,1])$ is expressed by means of a nested sequence of closed subspaces $\mathbf{V}_{0}\subset\mathbf{V}_{1}\subset\dots\subset\mathbf{V}_{j}\subset\dots\subset\mathbf{L}^2([0,1])$, where each subspace $\mathbf{V}_j$ corresponding to resolution $j$ is given by $\mathbf{V}_j = \mathrm{span}\{\phi_{j,k}(y): k = 0,\dots,2^{j}-1\}$.

Here, the functions

$$ \phi_{j,k}(y)=2^{j/2}\phi(2^{j}y-k) $$

are dilations and translations of a scaling function $\phi(y):[0,1]\rightarrow\mathbb{R}$, and $\phi(y)$ is such that the closure of the union of $\mathbf{V}_{j}$, $\overline{\bigcup_{k=0}^{\infty}\mathbf V_k}$, is dense in $\mathbf{L}^2([0,1])$.

Let the wavelet subspace $\mathbf W_j$ denote the orthogonal complement of $\mathbf{V}_{j}$ in $\mathbf{V}_{j+1}$, that is $\mathbf V_{j+1} = \mathbf{V}_{j} \oplus \mathbf W_j$ and $\mathbf{V}_{j}\bot \mathbf W_j$. It can be shown that $\mathbf W_j = \mathrm{span}\{\varphi_{j,k}(y):k = 0,\dots,2^{j}-1\}$ where $\varphi_{j,k}(y)$ is generated from dilation and translation of a so-called mother wavelet function $\varphi(y):[0,1]\rightarrow\mathbb{R}$,

$$ \varphi_{j,k}(y) = 2^{j/2}\varphi(2^{j}y-k). $$

By the construction of wavelet subspaces $\mathbf W_j$, it is straightforward to see that $\mathbf V_{j} = \mathbf V_0 \oplus\left(\bigoplus_{k=0}^{j-1}\mathbf W_k\right)$, and consequently $\mathbf V_0 \oplus\left(\bigoplus_{k=0}^{\infty}\mathbf W_k\right)=\mathbf{L}^2([0,1])$. Therefore, any function $u(y)\in \mathbf{L}^2([0,1])$ admits an orthogonal decomposition of the form

$$ u(y) = \tilde{\alpha}_{0,0}\phi_{0,0}(y) + \sum_{j=0}^{\infty}\sum_{k=0}^{2^{j}-1} \alpha_{j,k} \varphi_{j,k}(y), $$

where $\tilde{\alpha}_{0,0} = \langle u,\phi_{0,0} \rangle_{\mathbf{L}^2([0,1])}$, $\alpha_{j,k}= \langle u,\varphi_{j,k} \rangle_{\mathbf{L}^2([0,1])}$, and the inner-product $\langle u,v \rangle_{\mathbf{L}^2([0,1])} = \int_{0}^{1}u(y)\,v(y)\,dy$.

For the interest of notation, we rewrite the above expansion of $u(y)$ as

$$ u(y) = \sum_{i=1}^{\infty}\alpha_{i} \psi_{i}(y), $$

in which we establish a one-to-one correspondence between elements of the basis sets $\{\phi_{0,0},\varphi_{j,k}:k = 0,\dots,2^{j}-1, j = 0,\dots,\infty \}$ and $\{\psi_i:i = 1,\dots,\infty\}$.

Multiwavelet basis in one dimension

In the present study, we adopt the multiresolution basis of Alpert, where multiple scaling functions $\{\phi_i(y): i = 0,\dots,m-1\}$ are used to construct the polynomial space $\mathbf{V}^{m}_{0}$.

Specifically, we choose $\phi_i(y)$ as the Legendre polynomial of degree $i$ defined on $[0,1]$. An orthonormal piecewise polynomial basis $\{\varphi_i(y)\in\mathbf{L}^2([0,1]): i = 0,\dots,m-1\}$ of $\mathbf{W}^{m}_{0}$ is also established such that

$$ \int_{0}^{1}\varphi_i(y)y^l dy = 0,\qquad i,l=0,\cdots,m-1. $$

Let $\mathbf{U}_z = \{u(y)\in\mathbf{L}^2([-1,1]): \int^{1}_{-1}u(y)\,y^{z}\,dy=0\}$ represent the subspace of functions supported on the interval $[-1,1]$ with $z$ vanishing moments.

Given a maximum polynomial degree $m$, we construct $\{f_i:\,i=0,\dots,m-1\}$ in steps:

• Start from monomials of increasing degree, discontinuous and opposite in sign on the sub-intervals $[-1,0)$ and $[0,1]$.

• Require that $f_i\in\mathbf{U}_{j}$, $j=0,\dots,i+m-1$.

• Enforce orthogonality, i.e., $\langle f_i,f_j\rangle_{\mathbf{L}^2[(0,1)]}=\delta_{ij}$, $i,j = 0,\dots,m-1$, where $\delta_{ij}$ is the Kronecker delta.

• Finally compute $\{\varphi_i(y): i = 0,\dots,m-1\}$ by normalization $\{f_i:\,i=0,\dots,m-1\}$ on the interval $[0,1]$.

The multiwavelet basis functions $\varphi_{j,i,k}(y)$, hence the multiwavelet subspaces $\mathbf{W}^{m}_j$, are generated by dilations and translations of $\{\varphi_i(y): i = 0,\dots,m-1\}$, that is

$$ \varphi_{j,i,k}(y) = 2^{j/2}\varphi_i(2^j y-k),\quad i=0,\cdots,m-1,\; k=0,\cdots,2^j-1. $$

With certain additional constraints on $\varphi_i(y)$, the resulting basis is unique (up to the sign) and provides a generalization of the Legendre and Haar basis. In particular, the Legendre polynomials may be obtained by limiting the resolution to $j=0$, while Haar wavelets are obtained by setting $m=1$.

You can play with the maxOrder, maxDetailLevel, use random or equally spaced sampling to see how they affect the generated approximant.