Regression-linked estimation for a known graph
Section of Graphical models, part of Statistical learning
Regression-linked estimation for a known graph
| Type | Section |
| Taught in | Graphical models |
| Part | Statistical learning |
Several methods fit a Gaussian graphical model with a known partial covariance structure, for example the iterative proportional fitting procedure of Speed and Kiiveri (1986). The approach below can be extended to model selection using a form of the lasso.
The constrained likelihood
Main article: Constrained likelihood for a graph
For iid observations x_1,\ldots,x_n from \operatorname{N}_p(\mu, K^{-1}) the likelihood is
L=\prod_{i=1}^{n} \frac{|K|^{1/2}}{(2\pi)^{p/2}} \exp\Bigl(-\tfrac12(x_i-\mu)^\top K(x_i-\mu)\Bigr) =\frac{|K|^{n/2}}{(2\pi)^{np/2}} \exp\Biggl(-\tfrac12\sum_{i=1}^{n}(x_i-\mu)^\top K(x_i-\mu)\Biggr),
and the log-likelihood is
\ell=\frac{n}{2}\log|K|-\frac{np}{2}\log(2\pi) -\frac12\sum_{i=1}^{n}(x_i-\mu)^\top K(x_i-\mu).
Maximise it subject to constraints on K: edges missing from the graph correspond to zeros in K. If E is the edge set then k_{jk}=0 for j\neq k with (j,k)\notin E. The domain is K\succ 0. Diagonal entries are not constrained to zero. The second term is not a function of the parameters, and the likelihood is maximised for \mu at \hat\mu=\bar x regardless of K. Divide by n/2 without affecting the maximisation, so that one maximises \log|K|-\frac1n\sum_i(x_i-\bar x)^\top K(x_i-\bar x).
To simplify the later notation, let
S=\frac1n\sum_{i=1}^{n}(x_i-\bar x)(x_i-\bar x)^\top,
the covariance with divisor n. The unbiased matrix uses divisor n-1. In this section the symbol S means the 1/n version. Since a scalar equals its trace, and by the cyclic property of the trace,
\frac1n\sum_i(x_i-\bar x)^\top K(x_i-\bar x) =\operatorname{tr}\Biggl(\frac1n\sum_i(x_i-\bar x)(x_i-\bar x)^\top K\Biggr) =\operatorname{tr}(SK).
The problem is to maximise
g(K)=\log|K|-\operatorname{tr}(SK)
subject to the fixed zeros of K. With Lagrange multipliers \gamma_{jk} for the constraints,
g_C(K,\Gamma)=\log|K|-\operatorname{tr}(SK)-\sum_{j\neq k,\,(j,k)\notin E}\gamma_{jk} k_{jk}.
Here \Gamma is symmetric, zero on the diagonal and on allowed edges, and potentially nonzero on missing edges. A full-matrix trace inner product is used, so each off-diagonal pair appears twice. The symmetric multiplier convention absorbs that factor consistently. The derivatives with respect to symmetric K are K^{-1}=\Sigma for the log determinant and S for the trace term, giving \nabla g_C=\Sigma-S-\Gamma.
Using W for the estimate of \Sigma and \Theta for the estimate of K, so W\Theta=I, the estimating equations are
W-S-\Gamma=0.
Reading this entrywise: for every edge (j,k)\in E and for the diagonal, \gamma_{jk}=0 and so w_{jk}=s_{jk}. The maximum-likelihood estimate of the covariance reproduces the sample covariances on the edges and the diagonal exactly. For non-edges, w_{jk} is whatever value makes (W^{-1})_{jk}=0.
Solving by regression
Main article: Fitting a graph by regression
Isolate the last row and column of W and of its inverse \Theta. Any variable can be moved to the last position. Write A for the first p-1 indices and
\begin{pmatrix} W_{AA} & w_{Ap} \\ w_{Ap}^\top & w_{pp} \end{pmatrix} \begin{pmatrix} \Theta_{AA} & \theta_{Ap} \\ \theta_{Ap}^\top & \theta_{pp} \end{pmatrix} =\begin{pmatrix} I & 0 \\ 0^\top & 1 \end{pmatrix}.
The top p-1 entries of the last column of the identity give W_{AA}\theta_{Ap}+w_{Ap}\theta_{pp}=0. Since \theta_{pp} is a scalar,
w_{Ap}=-W_{AA}\theta_{Ap}/\theta_{pp}=W_{AA}\beta, \qquad \beta=-\theta_{Ap}/\theta_{pp}.
The top p-1 entries of the last column of the estimating equation are w_{Ap}-s_{Ap}-\gamma_{Ap}=0, which becomes
W_{AA}\beta-s_{Ap}-\gamma_{Ap}=0.
These p-1 equations are the estimating equations for a constrained regression of X_p on the other variables, with W_{AA}, the current estimate of \Sigma_{AA}, used in place of S_{AA}. The unconstrained normal equations of the regression of X_p on X_A are S_{AA}\beta=s_{Ap}.
On a complete sample those coefficients are exactly the ordinary least-squares coefficients, whenever every entry of S_{AA} and of s_{Ap} is computed with the same divisor. The factor n or n-1 cancels in the product S_{AA}^{-1}s_{Ap}. The adjective attached to the 1/n covariance is a statement about the covariance estimator. It does not move the fitted conditional-mean line relative to least squares.
The link with regression is the conditional normal. For a multivariate normal vector the conditional mean of X_p given X_A=x_A is \mu_p+\Sigma_{pA}\Sigma_{AA}^{-1}(x_A-\mu_A), the population multiple regression of X_p on X_A with \beta=\Sigma_{AA}^{-1}\Sigma_{Ap}. The block-inverse formula, with the last variable placed first, gives K_{Ap}=-k_{pp}\Sigma_{AA}^{-1}\Sigma_{Ap}, that is
\Sigma_{AA}^{-1}\Sigma_{Ap}=-K_{Ap}/k_{pp},
the same as \beta above except that the latter is an estimate. Estimating \Sigma or K can be viewed inside the linear regression framework, constrained by the known zeros of \Theta.
Suppose X_p has q neighbours. Its p-1-q missing edges force the same coordinates of \theta_{Ap}, and hence of \beta, to zero. On the q allowed edges the multipliers are zero, so those equations identify the q free regression coefficients. The remaining equations determine the missing-edge multipliers. Retain the neighbour rows and columns to obtain
W_{AA}^{*}\beta^{*}-s_{Ap}^{*}=0, \qquad \hat\beta^{*}=\bigl(W_{AA}^{*}\bigr)^{-1}s_{Ap}^{*},
where the starred blocks keep only the rows and columns of the neighbours of p. Pad \hat\beta^{*} with zeros to length p-1 to get \hat\beta. From the bottom-right entry of the block identity, w_{Ap}^\top\theta_{Ap}+w_{pp}\theta_{pp}=1, and with \theta_{Ap}=-\hat\beta\,\theta_{pp} one gets 1/\theta_{pp}=w_{pp}-w_{Ap}^\top\hat\beta, where w_{pp}=s_{pp} because the diagonal is treated as an edge.
Definition 1 (Maximum likelihood for a known graph)
- Initialise W=S. This presentation assumes S\succ 0.
- Repeat for j=1,\ldots,p until convergence:
- Partition W into part A, all but the jth row and column, and part B, the jth row and column.
- Solve W_{AA}^{*}\beta^{*}-s_{AB}^{*}=0 for the unconstrained edge parameters \beta^{*}. Obtain \hat\beta by padding \hat\beta^{*} with zeros in the positions of the missing edges.
- Update w_{AB}=W_{AA}\hat\beta and the corresponding row w_{BA}=w_{AB}^\top.
- In the final cycle for each j, solve for the estimated precision: \hat\theta_{AB}=-\hat\beta\,\hat\theta_{BB} with 1/\hat\theta_{BB}=s_{BB}-w_{AB}^\top\hat\beta.
The algorithm is the derivation applied to each variable in turn after moving it to the last position. Every update must maintain symmetry and positive definiteness, and convergence should be assessed through both parameter changes and the likelihood equations. A known graph can permit a maximum-likelihood estimate with singular S, but that case requires an appropriate feasible start and an existence check. The inverse operations above cannot simply be applied to singular blocks.
Example 1 (Fitting a chain) Question. Forty observations from the four-variable autoregressive model give the 1/n sample covariance
S=\begin{pmatrix} 0.8088 & 0.3166 & 0.2911 & 0.1228 \\ 0.3166 & 0.5490 & 0.4638 & 0.2318 \\ 0.2911 & 0.4638 & 0.9569 & 0.5472 \\ 0.1228 & 0.2318 & 0.5472 & 0.8223 \end{pmatrix}.
Fit the chain 1-2-3-4.
Solution. In the sweep for j=1 the only neighbour of variable 1 is 2, so the reduced regression of X_1 on X_2 has
\hat\beta^{*}=s_{12}/w_{22}=0.3166/0.5490=0.577,
and the update sets w_{12}=\hat\beta^{*} w_{22}=s_{12}, so the edge is reproduced exactly, while w_{13}=\hat\beta^{*} w_{23} and w_{14}=\hat\beta^{*} w_{24}. Each non-edge covariance is filled in from the covariances along the path through the neighbour. The sweeps for the other variables do the same. The constrained covariance, rounded to four decimals, is
W=\begin{pmatrix} 0.8088 & 0.3166 & 0.2675 & 0.1530 \\ 0.3166 & 0.5490 & 0.4638 & 0.2652 \\ 0.2675 & 0.4638 & 0.9569 & 0.5472 \\ 0.1530 & 0.2652 & 0.5472 & 0.8223 \end{pmatrix}, \qquad W^{-1}=\begin{pmatrix} 1.597 & -0.921 & 0 & 0 \\ -0.921 & 3.616 & -1.495 & 0 \\ 0 & -1.495 & 2.412 & -1.123 \\ 0 & 0 & -1.123 & 1.963 \end{pmatrix}.
The diagonal and the three edge entries of W equal those of S. For example w_{23}=s_{23}=0.4638. The non-edge entries are filled in so that the inverse has exact zeros where the graph has no edges: w_{13}=0.2675 against the sample value s_{13}=0.2911, w_{14}=0.1530 against s_{14}=0.1228, and w_{24}=0.2652 against s_{24}=0.2318. This is the Gaussian maximum-likelihood estimate for the constrained model. For example
w_{13}=w_{12}w_{23}/w_{22}=0.3166\times 0.4638/0.5490=0.2675,
the relation implied by X_1\perp X_3\mid X_2 in a chain, and similarly w_{14}=w_{12}w_{24}/w_{22}=0.1530 and w_{24}=w_{23}w_{34}/w_{33}=0.2652.
Used in
See also
Undirected graph · Hammersley–Clifford · A zero precision entry is conditional independence · Maximum likelihood for a known graph · Graphical lasso · Markov graphs · Graph terminology · Conditional independence