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

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?

• Have you followed the reference "Yserentant "? – faleichik Feb 4 '12 at 20:51
• Yes, but it only makes a statement about $S$, not about the decomposition. – rsp1984 Feb 4 '12 at 23:03
• Correction: "Yserentant" does actually make a statement about $S_i$ since he makes the presentation for $L=2$, implying that $S = S_1$. Will try to understand his approach. – rsp1984 Feb 4 '12 at 23:37

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  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  give a more detailed derivation of the structure (Eq. 9.29 - 9.32), however they use a rather unpractical inverse definition.

Combining both  and  I found the following to describe the decomposition structure pretty well: Following  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.

• Rafael, how important is the vector notation? It seems like using the \vec command doesn't necessarily render well (at least on my screen; perhaps it looks pretty on yours). If expressing the idea that it's a vector is important, could you suggest (or use) alternate notation? – Geoff Oxberry Feb 5 '12 at 19:59
• Happy to change notation. What other ways are there to express vectors here? – rsp1984 Feb 5 '12 at 20:04
• Some people follow the math community's convention and don't use any special notation to denote vectors. Others seem to follow the engineering convention of using bold (that is, \mathbf) to denote matrices and vectors. Feel free to use your preference. If the vector arrows rendered nicely, then I think your notation would be fine, but as it stands, some of the vector arrows are far to the left of where they should be, giving them kind of a disembodied look. – Geoff Oxberry Feb 5 '12 at 20:09