Exact finite Gaussian block elimination

The exact Gaussian block-elimination algorithm has input:

• an integer $N\ge 1$ and a real matrix $K_0\in {\mathbb{R}}^{N\times N}$ satisfying $K_0^{\mathsf{T}}=K_0$ and $v^{\mathsf{T}}{K}_{0}v>0$ for every nonzero $v\in {\mathbb{R}}^N$. In the finite Gaussian convention, $K_0$ is the precision because it appears in the density proportional to $\exp (-z^{\mathsf{T}}{K}_{0}z\mathrm{/}2)$;
• an integer $q\ge 0$ and nested retained coordinate sets $R_0=\{1,\dots ,N\}\supseteq R_1\supseteq \cdots \supseteq R_q,$ each given in increasing coordinate order.

It starts at stage $0$ with the ordered coordinates $R_0$ and matrix $K_0$. For each $j=1,\dots ,q$, let $E_j=R_{j-1}\setminus R_j$ be the coordinates to eliminate at that stage.

1. If $E_j=\mathrm{\varnothing }$, set $K_j=K_{j-1}$ and ${\mathrm{\ell }}_j=0$.
2. If $E_j\mathrm{\ne }\mathrm{\varnothing }$, reorder the current coordinates as the increasing list of $R_j$ followed by the increasing list of $E_j$, and write the correspondingly permuted matrix $K_{j-1}=\left(\begin{array}{cc}{A}_j& B_j\\ B_j^{\mathsf{T}}& C_j\end{array}\right)\mathrm{.}$ Solve the nonsingular system $C_j{X}_j=B_j^{\mathsf{T}}$, set $K_j:=A_j-B_j{X}_j,{\mathrm{\ell }}_j:=\frac{\mathrm{\mid }{E}_j\mathrm{\mid }}{2}\ln (2\pi )-\frac{1}{2}\ln (\det C_j)\mathrm{.}$ Interpret $K_j$ in the increasing coordinate order on $R_j$.

The output is $(K_q,{\mathrm{\ell }}_1,\dots ,{\mathrm{\ell }}_q)$, with the accumulated log normalization ${\sum }_{j=1}^q{\mathrm{\ell }}_j$. If $q=0$, this output means $K_0$ and an empty list of increments. When $R_j=\mathrm{\varnothing }$, $K_j$ is the unique empty $0\times 0$ matrix and its determinant is defined to be $1$. This is an exact-arithmetic procedure on finite matrices.

Finite Gaussian probability with a positive precision matrix

Positive Schur complements for finite Gaussian precisions

The natural logarithm as the inverse of the exponential function