Constructing Precision Matrices based on Correlated Gaussian Samples

Introduction

Here is a draft of a method to construct precision matrices based on correlated Gaussian samples with mean zero and variance one.

Method

Let $y_t$ be a random vector of size $J$ and $t\in \{1,\dots ,T\}$, where $T$ is the number of temporal replicates. We assume that $y_t\sim \mathcal{N}(0,Q^{-1})$ and that the marginal variance of the $i$-th element of $y_t$, $y_{i,t}$, is $1$. This means that the diagonal of $Q^{-1}$ is a vector of ones, and that

$$E(y_t)=0,\mathrm{cov}(y_t)=Q^{-1}.$$

Furthermore, it is assumed that Q is a sparse precision matrix. Using the properties of Gaussian conditional distributions, we have

$$E(y_{i,t}|y_{-i,t})=-Q_{i,i}^{-1}\sum _{j\in {\mathcal{A}}_i,j\ne i}{Q}_{i,j}{y}_{j,t},$$

$$\mathrm{Prec}(y_{i,t}|y_{-i,t})=Q_{i,i}=(\mathrm{var}({\mathrm{y}}_{\mathrm{i},\mathrm{t}}|{\mathrm{y}}_{-\mathrm{i},\mathrm{t}}){)}^{-1}={\tau }_i^{-2},$$

where ${\mathcal{A}}_i$ is the set containing the neighbors of site $i$, i.e. the sites that are such that $Q_{i,j}\ne 0$ if $j\in {\mathcal{A}}_i$.

Assume that we have realizations of $y_1,\dots y_t$ that can be used to infer the precision matrix $Q$. We set up a regression model to estimate the non-zero elements of $Q$. Here, we consider $y_{i,t}$ as a realization, i.e. as an observation. The regression model for each site, $i$, will be

$$y_{i,t}=\sum _{j\in \mathcal{A},j\ne i}{\beta }_{i,j}{y}_{j,t}+{\epsilon }_{i,t},t\in \{1,\dots ,T\}.$$

At each site $i$, we estimate the parameter vector ${\beta }_i$ with

$${\hat{\beta }}_i=(X_i^T{X}_i{)}^{-1}{X}_i^T{y}_i,$$

where

$$X_i=\left(\begin{array}{ccc}{y}_{j_{1,i},1}& \dots & y_{j_{m,i},1}\\ ⋮& ⋮& ⋮\\ y_{j_{1,i},T}& \dots & y_{j_{m,i},T}\end{array}\right),$$

and $y_{j_{l,i},1}$ is the $l$-th neighbor oy $y_{i,t}$ at time $t$. The variance of ${\epsilon }_{i,t}$ is ${\tau }_i^2$ and it is estimated with

$${\hat{\tau }}_i^2=T^{-1}(y_i-X_i{\hat{\beta }}_i{)}^T(y_i-X_i{\hat{\beta }}_i).$$

The next step is to transform ${\hat{\beta }}_i$ and ${\hat{\tau }}_i^2$ such that they give estimates of the elements of $Q$, namely

$${\hat{Q}}_{i,j}=\{\begin{array}{l}-{\hat{\tau }}_i^2{\hat{\beta }}_{i,j},\text{if }i\ne j,\\ {\hat{\tau }}_i^2,\text{ if }i=j,\end{array}$$

where ${\hat{\beta }}_{i,j}$ is the $j$-th element og ${\hat{\beta }}_i$. Let $\hat{B}$ be a matrix with $(i,j)$-th element ${\hat{\beta }}_{i,j}$. Note that ${\hat{\beta }}_{i,i}=0$, and thus ${\hat{B}}_{i,i}=0$. Furthermore let $\hat{K}$ be a diagonal matrix such that

$$\hat{K}=\mathrm{diag}({\hat{\tau }}_1^{-2},\dots ,{\hat{\tau }}_J^{-2}).$$

An estimate of Q can now be presented as

$$\hat{Q}=\hat{K}(I+\hat{B}),$$

where $I$ is an identity matrix of size $J$.

We have to make sure that $\hat{Q}$ is symmetric. This can be achieved by setting

$$\tilde{Q}i,j={\tilde{Q}}_{j,i}=\frac{1}{2}({\hat{\tau }}_i^{-2}{\hat{\beta }}_{i,j}+{\hat{\tau }}_j^{-2}{\hat{\beta }}_{j,i}),$$

and defining new regression parameters ${\tilde{\beta }}_{i,j}$ that are such that

$${\hat{\tau }}_i^{-2}{\tilde{\beta }}_{i,j}={\tilde{Q}}_{ij}={\tilde{Q}}_{j,i}={\hat{\tau }}_j^{-2}{\tilde{\beta }}_{j,i},$$

which gives

$${\tilde{\beta }}_{i,j}={\hat{\tau }}_i^2{\tilde{Q}}_{i,j},{\tilde{\beta }}_{j,i}={\hat{\tau }}_j^2{\tilde{Q}}_{i,j},$$

and let $\tilde{Q}$ and $\tilde{B}$ be the matrices containing the ${\tilde{Q}}_{i,j}$’s and the ${\tilde{\beta }}_{i,j}$’s.

We can not be sure of $\tilde{Q}$ being positive definite. One way to check whether the matrix is positive definite or not, is to compute the Cholesky decomposition of $\tilde{Q}$, that is, $\tilde{Q}=L{L}^T$, and check whether all the diagonal elements of L are positive. If the matrix $\tilde{Q}$ is invertible then it is more likely that it is positive definite, while if $\tilde{Q}$ is not invertible then it is not positive definite. The estimated precision matrix, $\tilde{Q}$, is invertible if $(i+\tilde{B})$ is invertible, where $\tilde{Q}=\hat{K}(I+\tilde{B})$. Strictly diagonally dominant matrices are invertible. In general, the $n\times n$ A, with elements ${\left\{a_{i,j}\right\}}_{i,j}$, is strictly diagonally dominant if

$$|a_{i,i}|>\sum _{j\ne i}|a_{i,j}|,1\le i\le n.$$

The matrix $(I+\tilde{B})$ is strictly diagonally dominant if

$$1>\sum _{j\in {\mathcal{A}}_i}|{\tilde{\beta }}_{i,j}|,1\le i\le J,$$

for all $i\in \{1,\dots ,J\}$. Alternatively, ${\lambda }_i\in (0,1)$ is found for each $i$ to tune $\tilde{Q}$ such taht it is strictly diagonally dominant, using

$${\hat{\tau }}_i^{-2}>{\lambda }_i\sum _{j\in {\mathcal{A}}_i,j\ne i}|{\tilde{Q}}_{i,j}|,1\le i\le J,$$