Correctness of finite Gaussian block elimination

Initially $K_0$ is positive definite. At any later stage with nonempty $E_j$, the current $K_{j-1}$ is positive definite by induction, so its principal block $C_j$ is positive definite. The Schur identity in [F2] then gives exactly the algorithm's $K_j=A_j-B_j{C}_j^{-1}{B}_j^{\mathsf{T}}$ and proves it positive definite when $R_j$ is nonempty. If $R_j$ is empty, the empty-matrix convention applies. Stages with $E_j=\mathrm{\varnothing }$ leave $K_j=K_{j-1}$. The algorithm is consequently well-defined at every one of its finitely many stages.

Fix one stage with $E_j\mathrm{\ne }\mathrm{\varnothing }$ and fix its retained coordinates $x$. By step 1.1 the current matrix is well-defined and positive definite. Reordering the coordinates is a linear permutation with Jacobian determinant of absolute value $1$, so the change-of-variables formula in [F3] preserves the integral. With eliminated coordinates denoted by $y$, [F2] gives ${\left(\begin{array}{c}x\\ y\end{array}\right)}^{\mathsf{T}}{K}_{j-1}\left(\begin{array}{c}x\\ y\end{array}\right)=x^{\mathsf{T}}{K}_{j}x+(y+C_j^{-1}{B}_j^{\mathsf{T}}x{)}^{\mathsf{T}}{C}_j(y+C_j^{-1}{B}_j^{\mathsf{T}}x)\mathrm{.}$ For fixed $x$, the translation $y\mapsto y+C_j^{-1}{B}_j^{\mathsf{T}}x$ is a $C^1$ diffeomorphism of Euclidean space with Jacobian $1$. The change-of-variables formula [F3] and the Gaussian normalizer [F1] therefore give ${\int }_{{\mathbb{R}}^{E_j}}{e}^{-\frac{1}{2}(x,y{)}^{\mathsf{T}}{K}_{j-1}(x,y)}dy=\frac{(2\pi {)}^{\mathrm{\mid }{E}_j\mathrm{\mid }\mathrm{/}2}}{\sqrt{\det C_j}}{e}^{-\frac{1}{2}{x}^{\mathsf{T}}{K}_{j}x}=e^{{\mathrm{\ell }}_j}{e}^{-\frac{1}{2}{x}^{\mathsf{T}}{K}_{j}x}\mathrm{.}$ The last equality follows from the definition of ${\mathrm{\ell }}_j$ and positivity of $\det C_j$. If $E_j=\mathrm{\varnothing }$, the same identity holds with factor $1=e^{{\mathrm{\ell }}_j}$ by definition.

Apply step 2.1 successively for the finite stages. At each stage the integrand is nonnegative, so Tonelli [F3] permits the iterated eliminations in the stated order. Multiplying the one-stage factors yields, for each fixed final retained vector $x$, ${\int }_{{\mathbb{R}}^{{\bigcup }_j{E}_j}}{e}^{-\frac{1}{2}\tilde{z}(x,y{)}^{\mathsf{T}}{K}_0\tilde{z}(x,y)}dy=\left({\prod }_{j=1}^q{e}^{{\mathrm{\ell }}_j}\right)e^{-\frac{1}{2}{x}^{\mathsf{T}}{K}_{q}x}=e^{{\sum }_j{\mathrm{\ell }}_j}{e}^{-\frac{1}{2}{x}^{\mathsf{T}}{K}_{q}x},$ using exponential addition in [F4]. If ${\bigcup }_j{E}_j=\mathrm{\varnothing }$, including $q=0$, the integral is with respect to the unit mass on ${\mathbb{R}}^0$ and the equality is the identity with an empty product.

If no coordinates are eliminated, then every $E_j=\mathrm{\varnothing }$, $K_q=K_0$, and ${\mathrm{\ell }}_j=0$, so the normalizer identity is immediate. If $R_q=\mathrm{\varnothing }$, then ${\bigcup }_j{E}_j=R_0$ has $N\ge 1$ coordinates; evaluating step 3.1 at the unique empty retained vector gives $Z_{K_0}=e^{{\sum }_j{\mathrm{\ell }}_j}$, while $Z_{K_q}=1$. In the remaining case both $R_q$ and ${\bigcup }_j{E}_j$ are nonempty, so Tonelli [F3] and the Euclidean product-measure identity in [F3] allow integration of step 3.1 over the retained vector and give $Z_{K_0}=e^{{\sum }_j{\mathrm{\ell }}_j}{Z}_{K_q}$. The nonempty $K_q$ is positive definite by step 1.1, so [F1] makes $Z_{K_q}$ finite and positive. In each case, dividing step 3.1 by $Z_{K_0}$ gives the retained density $Z_{K_q}^{-1}{e}^{-x^{\mathsf{T}}{K}_{q}x\mathrm{/}2}$ when $R_q\mathrm{\ne }\mathrm{\varnothing }$, exactly ${\gamma }_{K_q}$; for empty $R_q$ the marginal is the unit mass ${\gamma }_{\mathrm{\varnothing }}$. Taking natural logarithms of the positive normalizer identity and using [F4] proves the logarithmic formula, including the empty-space case.