Proof of Mean-Square Limits of Gaussian Random Vectors are Gaussian

Each component $X^{k}_{i}$ is a Gaussian random variable, since selecting a coordinate of a Gaussian random vector is an affine transformation (Affine Transformations of Gaussian Random Vectors are Gaussian); in particular each $X^k_i$ is square-integrable by Square-Integrability, Moments, and Covariance Matrix of a Gaussian Random Vector.

Step 1 (Componentwise limits and moment convergence). By Mean-Square Limits of Gaussian Random Variables are Gaussian applied to each coordinate, every $X_i$ is Gaussian with $\mu^{k}_{i}:=\mathbb{E}[X^k_i]\to\mathbb{E}[X_i]=:\mu_i$ and $\operatorname{Var}(X^k_i)\to\operatorname{Var}(X_i)$. Write $Y_i=X_i-\mu_i$ and $Y^{k}_{i}=X^{k}_{i}-\mu^{k}_{i}$; then

by the triangle inequality of Cauchy-Schwarz and Triangle Inequalities for the Mean-Square Norm. For the covariances, $\operatorname{Cov}(X^k_i,X^k_j)=\langle Y^k_i,Y^k_j\rangle_{2}$ and $\operatorname{Cov}(X_i,X_j)=\langle Y_i,Y_j\rangle_{2}$ with the mean-square inner product of Square-Integrable Random Variables and the Mean-Square Inner Product, and by the Cauchy-Schwarz inequality,

since $\lVert Y^k_j\rVert_2$ is bounded (a convergent sequence of norms). This proves the moment-convergence claims; it remains to show $(X_1,\dots,X_d)$ is a Gaussian random vector, and for this the final relabeling clause of the statement lets us discard finitely many initial indices whenever convenient.

Step 2 (Gram matrices). Define real $d\times d$ matrices $G$ and $G^{(k)}$ by their entries $G_{ij}=\langle Y_i,Y_j\rangle_2$ and $G^{(k)}_{ij}=\langle Y^k_i,Y^k_j\rangle_2$; they are symmetric by symmetry of the mean-square inner product (Square-Integrable Random Variables and the Mean-Square Inner Product), and $G^{(k)}\to G$ entrywise by Step 1. For every $c\in\mathbb{R}^{d}$, bilinearity of the mean-square inner product gives, with the dot product and matrix-vector product,

so all these matrices are positive semidefinite.

Step 3 (Selection of a positive definite block). Choose $I\subseteq\{1,\dots,d\}$ of maximal cardinality such that the submatrix $G_I=(G_{ij})_{i,j\in I}$ is positive definite (the empty set qualifies vacuously).

If $I=\emptyset$, then for each $i$ the singleton $\{i\}$ fails positive definiteness, i.e., $G_{ii}=\lVert Y_i\rVert_2^{2}=0$, so $P(X_i=\mu_i)=1$ by the null-equivalence clause of Square-Integrable Random Variables and the Mean-Square Inner Product; then $P\bigl(X_i=\mu_i+\text{(empty sum)}\bigr)=1$ for every $i$, which is a Gaussian representation with $m=0$ in the sense of Gaussian Random Vectors and Jointly Gaussian Random Variables, and we are done. So assume $I=\{i_1<\dots<i_r\}$ with $r\ge1$.

For $i\notin I$ we claim there are reals $(c^{(i)}_j)_{j\in I}$ with $Y_i=\sum_{j\in I}c^{(i)}_jY_j$ almost surely. Indeed, by maximality the submatrix over $I'=I\cup\{i\}$ is not positive definite; since it is positive semidefinite (its quadratic form is the restriction of that of $G$ to vectors supported on $I'$), there is a nonzero vector $c$ indexed by $I'$ with $c\cdot(G_{I'}c)=0$, i.e., $\lVert\sum_{j\in I'}c_jY_j\rVert_2^{2}=0$, whence $\sum_{j\in I'}c_jY_j=0$ almost surely. If $c_i=0$, then $c$ restricted to $I$ would be a nonzero vector annihilating the quadratic form of $G_I$, contradicting positive definiteness; so $c_i\ne0$ and $Y_i=\sum_{j\in I}(-c_j/c_i)Y_j$ almost surely.

Step 4 (Orthonormalization and independent standard normals). The $r\times r$ matrices $G^{(k)}_I$ are symmetric positive semidefinite and converge entrywise to the positive definite $G_I$. By Triangular Orthonormalization of a Positive Definite Gram Matrix there are $K$ and lower triangular matrices $T^{(k)}$ ($k\ge K$) and $T$, all with positive diagonals, with

and $T$ invertible with lower triangular inverse $S=T^{-1}$. Discarding the indices $k<K$ and relabeling (as allowed by the final clause of the statement), we may assume the $T^{(k)}$ are defined for all $k\in\mathbb{N}$.

For $1\le l\le r$ define

$(U^k_1,\dots,U^k_r)$ is a Gaussian random vector, being an affine image of $(X^k_1,\dots,X^k_d)$ (Affine Transformations of Gaussian Random Vectors are Gaussian). By linearity, $\mathbb{E}[U^k_l]=0$, and bilinearity of the mean-square inner product gives, for component indices $l,l'$,

(the Kronecker symbol: $1$ if $l=l'$ and $0$ otherwise). Thus the components of $(U^k_1,\dots,U^k_r)$ are jointly Gaussian, pairwise uncorrelated (Covariance of Square-Integrable Random Variables; the means vanish), each with variance $1$; by Pairwise Uncorrelated Jointly Gaussian Random Variables are Independent they are independent, and by clause 2 of Standardization and Cumulative Distribution Function of a Gaussian Random Variable (standardization with $\mu=0$, $\sigma=1$) each $U^k_l$ is standard normal.

Moreover $U^k_l\to U_l$ in mean square: by the triangle inequality of Cauchy-Schwarz and Triangle Inequalities for the Mean-Square Norm,

using $T^{(k)}\to T$ and boundedness of $\lVert Y^k_{i_m}\rVert_2$. By Markov's inequality applied to the nonnegative random variables $(U^k_l-U_l)^{2}$, mean-square convergence implies convergence in probability, so by Independence is Preserved by Limits in Probability the limits $U_1,\dots,U_r$ are independent; and by Mean-Square Limits of Gaussian Random Variables are Gaussian each $U_l$ is Gaussian with mean $\lim 0=0$ and variance $\lim 1=1$, hence standard normal by clause 2 of Standardization and Cumulative Distribution Function of a Gaussian Random Variable.

Step 5 (Assembling the representation). Since $ST=I_r$ pointwise as matrices, for each $m\in\{1,\dots,r\}$ the finite linear combinations satisfy, at every sample point,

Hence for $i=i_m\in I$: $X_i=\mu_i+\sum_{l=1}^{r}a_{il}U_l$ everywhere on $\Omega$, with $a_{il}=S_{ml}$. For $i\notin I$: on the almost sure event of Step 3, $X_i=\mu_i+\sum_{j\in I}c^{(i)}_jY_j=\mu_i+\sum_{l=1}^{r}a_{il}U_l$ with $a_{il}=\sum_{m}c^{(i)}_{i_m}S_{ml}$. The intersection of these finitely many almost sure events is almost sure (finite subadditivity of the measure $P$ applied to the complements). Therefore

with $U_1,\dots,U_r$ independent standard normal random variables. By Gaussian Random Vectors and Jointly Gaussian Random Variables, $(X_1,\dots,X_d)$ is a Gaussian random vector. $\square$