Conjugate Gradient with Hierarchical Basis Functions: How can the hierarchical base be decomposed?

Question

I'm trying to implement a Conjugate Gradient solver using Hierarchical Basis Functions, following this paper.

In section 3 the paper says that the hierarchical basis matrix $S$ can be decomposed into a "series of very sparse matrices $S = S_1S_2...S_{L-1}$ " but it doesn't tell the exact structure of these matrices $S_i$, nor in which way they are sparse.

Therefore: What is the structure of $S_i$ (I assume they are all similar in a way)? Why is $S_i$ sparse? What is the interpretation of applying $S_i$ to a vector?

Answer 1

Some research into the problem led me to two publications that do a fair job explaining the structure:

1. Bank, Dupont, and Yserentant (1987)
2. A book chapter by Pinksy, Malhotra, and Thompson (1994) (Page 183).

Whereas [1] gives a fairly concise matrix structure of $\mathbf{S}_1$ (Eq. 5.3, two hierarchy levels only), it lacks a more detailed explanation.
The authors of [2] give a more detailed derivation of the structure (Eq. 9.29 - 9.32), however they use a rather unpractical inverse definition.

Combining both [1] and [2] I found the following to describe the decomposition structure pretty well: Following [1] and initially taking a 2-Level approach, we can express $\mathbf{S}=\mathbf{S}_1$ as

$\mathbf{S} = \begin{pmatrix} I & 0 \\ \mathbf{R} & I \end{pmatrix}$

so that if we compute $\mathbf{x} = \mathbf{S}\mathbf{y}$, we get a representation of the function in terms of nodal basis functions, coming from the hierarchical representation $\mathbf{y} = \begin{pmatrix} \mathbf{y}_2 & \mathbf{y}_1 \end{pmatrix}^T$ where $\mathbf{y_2}$ are the coarse-level coefficients and $\mathbf{y}_1$ are the fine-level coefficients .

Since the hierarchical basis functions have larger influence regions on the coarse level, we can see that we need to incorporate this influence at the nodal level, which is what $\begin{pmatrix} \mathbf{R} & I \end{pmatrix}$ does.

$\mathbf{R}$ is itself a sparse matrix and its rows contain the values of the coarse-level basis functions at the nodal-variable-location - which is mostly 0, except for the locations where the influence of a coarse-level function overlaps with the nodal-variable-location.
A typical row looks like $\mathbf{R}_i = \begin{pmatrix} \mathbf{0} & a & b & \mathbf{0}\end{pmatrix}$ for a 1-dimensional problem.

In order to arrive at the final decomposition $\mathbf{S} = \mathbf{S_1S_2...S_{L-1}}$ we need to apply the above recursively:

A complete set of hierarchical coefficients looks like $\mathbf{y} = \begin{pmatrix} \mathbf{y_L} & \mathbf{y_{L-1}} & ... & \mathbf{y_1}\end{pmatrix}^T$.
We first get rid of the level-L-coefficients by expressing this part of the signal in terms of level L-1 -coefficients:

$\mathbf{y}' = \begin{pmatrix} I & 0 & 0\\ \mathbf{R} & I & 0 \\ 0 & 0 & I\end{pmatrix}\cdot\mathbf{y} = \mathbf{S}_{L-1}\cdot\mathbf{y}$

Now the coarsest basis is level L-1.
We can now, again, express this in a finer basis (L-2) using $\mathbf{y}'' = \mathbf{S}_{L-2}\cdot\mathbf{y}'$ and so on.
$\mathbf{S}_{L-2}$ looks exactly like $\mathbf{S}_{L-1}$ with the difference that $\mathbf{R}$ is a bit larger since there are more L-2 coefficients than L-1 coefficients.

In the end we will arrive at $\mathbf{y}^{L-1} = \mathbf{x} = \mathbf{S}_1\mathbf{y}^{L-2}$.
$\mathbf{S_1}$ will then look exactly like described above.