Finite-support realization of a covariance matrix ($E=LS$)

Every real symmetric matrix $\Sigma$ satisfying $a^T\Sigma a\ge 0$ for every real vector $a$ is the covariance matrix of a random vector taking only finitely many values. The stated inequality means that $\Sigma$ is a positive semidefinite matrix. One does not need a multivariate normal distribution to realize a specified second-moment geometry. Conversely, a proposed covariance matrix must be positive semidefinite, since $a^T\Sigma a$ is the variance of the projection onto any real vector $a$.

Here is an explicit construction. If $r$ is the matrix rank of $\Sigma$, orthogonal diagonalization gives a real $d\times r$ matrix $L$ withLet $S_1,\dots ,S_r$ be independent Rademacher random variables, each equal to $-1$ or $1$ with probability $1/2$, and put $S=(S_1,\dots ,S_r{)}^T$. Then $\mathbb{E}S=0$ and $\mathbb{E}[S{S}^T]=I_r$, soThe random vector $E$ has at most $2^r$ possible values. Adding any constant vector $\mu$ prescribes its mean without changing its covariance matrix. If $r=0$, the constant vector alone suffices.

For an original two-dimensional example chooseThe four equally likely values of $E$ areTheir coordinate means vanish. Their squared first coordinates average to $4$, their squared second coordinates to $(9+1+1+9)/4=5$, and their coordinate products to $(6-2-2+6)/4=2$. The displayed covariance matrix is therefore realized exactly, without sampling error or an asymptotic argument.

These four values are a complete probability model with masses $1/4$, not automatically an independent random sample whose population is unknown. Applying the usual divisor-three sample covariance matrix to the same four rows would give $4\Sigma /3$. NIST explicitly distinguishes that unbiased sample formula with divisor $n-1$; probability-weighted population moments use the assigned masses instead.^[ref]

The sign construction need not minimize the number of points. If a distribution has $m$ positive-mass support points, the centered points obey one nontrivial weighted relation, so their span has dimension at most $m-1$. Its covariance matrix consequently has matrix rank at most $m-1$, requiring $m\ge r+1$ when $r>0$. Matching a covariance matrix also does not match the full probability distribution: this finite model and a continuous Gaussian model can have identical means and covariances while assigning very different probabilities to nonlinear events.