Some research into the problem led me to two publications that do a fair job explaining the structure:
- Bank, Dupont, and Yserentant (1987)
- 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.
