Gaussian Vector – Understanding Rank-Deficient Transformation of Gaussian Vector

multivariate normal distributionnormal distribution

Suppose we have a gaussian $n$-dimensional vector $X \sim \mathcal{N}\left(\mu, \Sigma\right)$

I know that for a non-singular linear map $A$, its image $Y = AX$ will also be a gaussian vector with distribution $\mathcal{N}\left(A\mu, A\Sigma A^T\right)$.

I am concerned about the case of rank-deficient linear map $A$ (i.e. the matrix $A$ doesn't have full rank).

Will the $AX$ still be gaussian?

Best Answer

Suppose $A : \mathbb R^n \to \mathbb R^p$. $AX$ will be Gaussian, and we can use characteristic functions to see this.

$$ \varphi_Y(t) = \text E[e^{it^TY}] = \text E[e^{it^T(AX)}] = \text E[e^{i(A^Tt)^TX}] = \varphi_X(A^Tt) \\ = \exp\left(i \mu^TA^Tt - \frac 12 t^T A\Sigma A^Tt\right) $$ so $Y \sim \mathcal N(A\mu, A\Sigma A^T)$ as hoped. The one complexity here is that because $A\Sigma A^T$ is low rank, $Y$ does not have a density in the usual sense on $\mathbb R^p$. The column space of $A$ is a low dimension subspace of $\mathbb R^p$, and since $Y = AX$ we know that this subspace is the support of $Y$. Such a subspace has a measure of zero w.r.t. the Lebesgue measure on $\mathbb R^p$ so the requirement of absolute continuity fails, which is necessary for the Radon-Nikodym theorem.

Update: we can recover the rank of $A$ with probability one if we can take multiple samples. Let $\newcommand{\X}{\mathbf X}$$\X=[X_1\mid \dots \mid X_n]$ contain $n$ iid samples from $\mathcal N_n(\mu,\Sigma)$ and let $\newcommand{\Y}{\mathbf Y}$$\Y=[Y_1 \mid \dots \mid Y_n]$ where $Y_i = AX_i$ so $\Y = A\X$. The key insight is that $\X$ will be full rank almost surely, so $\text{rank}(\Y) = \text{rank}(A\X) = \text{rank}(A)$.

Here's an example in R:

set.seed(123)
n <- 15; p <- 23
r <- 12  # rank of A
mu <- rnorm(n)
# getting a random positive definite covariance matrix
Sigma <- matrix(rnorm(n*n),n,n); Sigma <- Sigma %*% t(Sigma)
X <- t(MASS::mvrnorm(n, mu, Sigma))
stopifnot(Matrix::rankMatrix(X) == n)
A <- matrix(rnorm(p*r), p, r) %*% matrix(rnorm(r*n), r, n)
stopifnot(Matrix::rankMatrix(A) == r)
Y <- A %*% X
stopifnot(Matrix::rankMatrix(Y) == r)
Related Question