I have recently started to delve into the construction and parameterization of orthogonal polynomial bases (e.g., for polynomial chaos expansion) and there are some things that aren't yet clear to me, or that I am not confident in understanding fully. I was thus hoping that you might confirm or clarify some uncertain points for me. In advance: sorry, this is a long post, but I hope most of it might be trivial for experienced folks, and I want to be thorough.
From what I have gathered so far, using a polynomial basis is similar to a Fourier expansion, just that instead of using sin functions, one uses polynomials of different degrees. So if I have a univariate variable $x \in \mathbb{R}$ and a function $f(x) \rightarrow \mathbb{R}$, I could for example approximate it with a polynomial of third degree $D=3$:
$$f(x) \approx a_3x^3 + a_2x^2+a_1x+a_0$$
which one could reformulate as follows:
$$f(x) \approx \begin{bmatrix} a_0 & a_1 & a_2 & a_3 \end{bmatrix} \begin{bmatrix} 1 \\ x \\ x^2 \\ x^3 \end{bmatrix} $$
where the first vector are the polynomial coefficients $a$ and the second vector is the polynomial basis, which one could interpret as a four-dimensional space in which each point corresponds to a specific set of four coefficients. One could also use a different polynomial basis, e.g., Hermite polynomials instead:
$$f(x) \approx \begin{bmatrix} a_0 & a_1 & a_2 & a_3 \end{bmatrix} \begin{bmatrix} 1 \\ x \\ x^2 -1 \\ x^3 - 3x \end{bmatrix} $$
which would of course affect the coefficients for this approximation as well. Is this correct so far?
Now if I extend this problem to more dimensions, say $N=2$, so that $X \in \mathbb{R}^N$ and $f(X) \rightarrow \mathbb{R}^N$ instead, and if I assume the dimensions are independent, I get (with a Hermite polynomial basis):
$$f(X) = \begin{bmatrix} f_1(x_1) \\ f_2(x_2) \end{bmatrix} \approx \begin{bmatrix} a_{1,3}(x_1^3-3x_1) + a_{1,2}(x_1^2-1) + a_{1,1}x_1 + a_{1,0} \\ a_{2,3}(x_2^3-3x_2) + a_{2,2}(x_2^2-1) + a_{2,1}x_2 + a_{2,0} \end{bmatrix}$$
which one could reformulate like this:
$$f(X) \approx \begin{bmatrix} a_{1,0} & a_{1,1} & a_{1,2} & a_{1,3} & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 &a_{2,0} & a_{2,1} & a_{2,2} & a_{2,3} \end{bmatrix} \begin{bmatrix} 1 \\ x_1 \\ x_1^2 -1 \\ x_1^3 -3x_1\\ 1 \\ x_2 \\ x_2^2 -1 \\ x_2^3 -3x_2\\ \end{bmatrix} $$
where I suppose one could combine the two zeroth-order dimensions of the polynomial basis:
$$f(X) \approx \begin{bmatrix} a_{1,0} & a_{1,1} & a_{1,2} & a_{1,3} & 0 & 0 & 0\\ a_{2,0} & 0 & 0 & 0 & a_{2,1} & a_{2,2} & a_{2,3} \end{bmatrix} \begin{bmatrix} 1 \\ x_1 \\ x_1^2 -1 \\ x_1^3 -3x_1\\ x_2 \\ x_2^2 -1 \\ x_2^3 -3x_2\\ \end{bmatrix} $$
Once again, is this correct so far?
Now if I assume the dimensions are not independent, things get more foggy for me. From what I think I understand is that in case of dependence between two dimensions, we expand the polynomial basis with the products of all dimension pairs, so we obtain (for brevity, monomial basis again, with $N=2$ and $D=2$:
$$basis = \begin{bmatrix} 1 \\ x_1 \\ x_2 \\ x_1x_2\\ x_1^2 \\ x_2^2 \\ x_1^2x_2 \\ x_2^2x_1 \\ x_1^2x_2^2 \\ \end{bmatrix} $$
Where each dimension in canonical space $n=1,2$ only has coefficients $a_{n,...}$ in dimensions of the basis which it is included in and is zero everywhere else (i.e., dimension $n=1$ has non-zero coefficients only in dimensions which are either $1$ or include $x_1$ in some form [particularly relevant for $N > 2$]). I have also come across 'truncations of the expansion' which have a 'lower triangular structure', which seems to remove the dimensions whose combined order is higher than a threshold, e.g.:
$$basis = \begin{bmatrix} 1 \\ x_1 \\ x_2 \\ x_1x_2\\ x_1^2 \\ x_2^2 \end{bmatrix} $$
where we only retained dimensions of summed $order \leq 2$, so that $x_1^2x_2$ ($order = 3$), $x_2^2x_1$ ($order = 3$), and $x_1^2x_2^2$ ($order = 4$) are omitted.
I assume this is what this 'lower triangular structure' is about? Where the left figure would correspond to the $N=2,D=3,independent$ case from above (omitting all cross-terms), and the right figure would correspond to the $N=2, D=2, dependent$ case discussed just now (blue marks terms which exist, white marks terms we omitted)?
Is this correct? Is there something I have missed?

The question is a little long and hard to follow (it would be better if it was posted as multiple different questions), but from what I can see this is what you want:
Yes, when you have a multi-variate system you need to multiply the polynomial bases together as you have described. Yes, this can create a really large basis is in the end. This phenomenon is known as the curse of dimensionality. If you want a multivariate polynomial basis of dimension = 2, you would be required to multiple your two families of polynomials together, and ensure that the indices of the polynomial terms add to being $\leq 2$
This can obviously explode in size very quickly, so in order to get around this we can consider structures like the one you propose, the lower triangular structure. This is all about trying to find the most efficient representaion using this multi-variate polynomial basis. Yes, if the two dimensions are independent, you wouldn't need to multiply the basis together. BUT, if they are not you will need to, and often not all the terms will contribute equally as well.
Finding this most efficient representation is no easy task, and is an active area of research. But it would seem like often a lot of higher order terms can be ignored with negligible modelling error. However, this is not a mathematical guarantee.
Hope this helps.