Here we fix the factor loading matrix \(\mathbf {B}\in \mathbb {R}^{n \times r_{\text{true}}}\) .
The latent factor \(\mathbf {f}\in \mathbb {R}^{r_{\text{true}}}\)  is drawn from a coordinate-independent, zero-mean and unit-variant Normal-inverse Gaussian (NIG) distribution [1]}, [2]}. This distribution is skewed and heavy-tailed, which is aligned with our interest mentioned in sec: introduction.
The Gaussian noise \(\mathcal {N}(\mathbf {0}, \sigma ^2\mathbf {I}_n)\)  is referred to \({\varepsilon } \in \mathbb {R}^n\) .
Note that the \(d\) -th moment of \(\mathbf {x}\)  roughly admits a low-rank Tucker decomposition format:
\(\begin{array}{lll}{M}^d\!\!&\approx \mathbb {E}\left(\mathbf {x}^{\otimes d}\right)\!\!&\approx \,\,\, \displaystyle \mathbb {E}\left((\mathbf {B}\mathbf {f})^{\otimes d}\right) \\&& \displaystyle = \,\,\, \mathbb {E}\left((\mathbf {R}^*\mathbf {f})^{\otimes d}\right)\cdot \left(\mathbf {Q}^{*\top }, \dots , \mathbf {Q}^{*\top }\right).\end{array}\) 
