# 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 \left\{1, \dots, T\right\}$$, 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, \qquad \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\neq i} Q_{i,j}y_{j,t},$

$\mathrm{Prec}(y_{i,t}|y_{-i, t}) = Q_{i,i} = (\mathrm{var(y_{i,t}|y_{-i, 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} \neq 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\neq i} \beta_{i,j}y_{j,t} + \varepsilon_{i, t}, \quad t\in \left\{1, \dots, T\right\}.$

At each site $$i$$, we estimate the parameter vector $$\beta_i$$ with

$\hat\beta_i = (X_i^TX_i)^{-1}X_i^Ty_i,$

where

$X_i = \begin{pmatrix} y_{j_{1, i}, 1} & \dots & y_{j_{m, i}, 1} \\ \vdots & \vdots & \vdots \\ y_{j_{1, i}, T} & \dots & y_{j_{m, i}, T} \end{pmatrix},$

and $$y_{j_{l, i}, 1}$$ is the $$l$$-th neighbor oy $$y_{i, t}$$ at time $$t$$. The variance of $$\varepsilon_{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{cases} -\hat\tau_i^2\hat\beta_{i, j}, \quad \text{if } i \neq j, \\ \hat\tau_i^2, \qquad \quad \text{ if } i = j, \end{cases}$

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}\left(\hat\tau_1^{-2}, \dots, \hat\tau_J^{-2}\right).$

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} = \frac12(\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}, \quad \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 = LL^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

$\vert a_{i, i}\vert > \sum_{j\neq i} \vert a_{i, j}\vert, \qquad 1\leq i \leq n.$

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

$1 > \sum_{j \in \mathcal A_i} \vert \tilde \beta_{i, j}|, \qquad 1 \leq i \leq J,$

for all $$i \in \left\{1, \dots, J \right\}$$. 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\neq i} \vert \tilde Q_{i, j} \vert, \qquad 1\leq i \leq J,$