Normally one would describe a Gaussian distribution in some dimension $n$ via some mean $\mu\in\mathbb{R}^n$ and its covariance matrix $\Sigma\in\mathbb{R}^{n\times n}$. I want to reduce the amount of variables needed, to describe this same exact function.
An easy first step is, to realize, that $\Sigma$ must be symmetric, reducing the number of variables to $n + \frac{n(n+1)}{2}$, however thats really the only nice reduction i can find.
Onto another approach. In most applications, covariances are often depicted as ovals around the mean, with some sheared axis. This is, where this next approach is motivated from. As can be seen here, one can view sheared ovals as just rotated ovals, with "orthogonal axis". Now if we think vaguely in terms of covariances, one might think, that we can find $n$ axis, that are orthogonal to one another, that describe the covariance. If this were true, we could describe a Gaussian distribution via
- some mean $\mu\in\mathbb{R}^n$
- some $n$ variances ($\in\mathbb{R}$) along major axis
- some rotation in $\mathbb{R}^n$, i.e. $n$ angles
which would reduce the number of variables needed to describe our distribution drastically from $n + \frac{n(n+1)}{2}$ to $3n$. So now to the question.
Is this possible, i.e. can we describe any $n$-dimensional Gaussian distribution via $n$ linearly independent vectors, and a variance along these vectors (and the mean $\mu$)?
The main remark is that it is not a rotation matrix that will transform the reference (orthogonal) axes into the proper (non orthogonal in general) axes of your ellipsoid. It is a shear matrix (you have used the term in your question): which has the form $I+N$ where $N$ is a strictly upper triangular matrix
$$\begin{pmatrix}1&a_{12}&a_{13}&...&a_{1n}\\ 0&1&a_{23}&...&a_{2n}\\ 0&\cdots&\cdots&\cdots&\cdots\\ 0&\cdots&\cdots&\cdots&\cdots\\ 0&0&0&1&a_{n-1,n}\\ 0&0&0&0&1\cdots \end{pmatrix}$$
(think to the orthonormalization process which uses such matrices at least implicitly).
Such a matrix uses $\dfrac{n(n-1)}{2}$ parameters.
Consequence : your count can be made exact :
The sum $2n+\dfrac{n(n-1)}{2}$ coincides indeed with the first count : $n+\dfrac{n(n+1)}{2}$.