Finite-support realization of a covariance matrix (source code)

= Finite-support realization of a covariance matrix
{title2=$E=LS$}

Every real <symmetric matrix> $\Sigma$ satisfying $a^T\Sigma a\geq0$ 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$ with
$$
\Sigma=LL^T.
$$
Let $S_1,\ldots,S_r$ be independent <Rademacher random variables>, each equal to $-1$ or $1$ with probability $1/2$, and put $S=(S_1,\ldots,S_r)^T$. Then $\mathbb ES=0$ and $\mathbb E[SS^T]=I_r$, so
$$
\boxed{E=LS,\qquad \mathbb EE=0,\qquad
\operatorname{Cov}(E)=L I_rL^T=\Sigma.}
$$
The <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 choose
$$
L=\begin{pmatrix}2&0\\1&2\end{pmatrix},\qquad
\Sigma=\begin{pmatrix}4&2\\2&5\end{pmatrix}.
$$
The four equally likely values of $E$ are
$$
(2,3),\quad (2,-1),\quad (-2,1),\quad (-2,-3).
$$
Their 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.<NIST and SEMATECH e-Handbook of Statistical Methods>[Section6.5.4.1, Mean Vector and Covariance Matrix, unbiased sample formula]{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\geq 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.