Undirected graphical models and sparse precision estimation

The state space $\{0,1{\}}^m\times \{0,1{\}}^n$ has $2^{m+n}$ elements. Every term $\exp (-E(v,h))$ is finite and positive, so their finite sum $Z$ is finite and positive; therefore $p$ is a probability mass function.

For fixed $v$, put $c_j=b_j+(W^{T}v{)}_j$. Then $\exp (-E(v,h))=\exp (a^{T}v){\prod }_j\exp (h_j{c}_j)$, and summing over binary $h$ gives $\exp (a^{T}v){\prod }_j(1+\exp (c_j))$. Thus $p(h\mid v)$ is the product of Bernoulli masses with success probabilities $\exp (c_j)\mathrm{/}(1+\exp (c_j))=\sigma (c_j)$.

For fixed $h$, put $d_i=a_i+(Wh{)}_i$. Then $\exp (-E(v,h))=\exp (b^{T}h){\prod }_i\exp (v_i{d}_i)$, and summing over binary $v$ gives $\exp (b^{T}h){\prod }_i(1+\exp (d_i))$. Hence $p(v\mid h)$ is the product of Bernoulli masses with success probabilities $\exp (d_i)\mathrm{/}(1+\exp (d_i))=\sigma (d_i)$.

For fixed observed $v$, $p(v)=Z^{-1}\exp (a^{T}v){\sum }_h\exp (b^{T}h+v^{T}Wh)$. Differentiating these finite sums gives ${\mathrm{\partial }}_{W_{ij}}\log p(v)=v_i{\mathbb{E}}_{p(h\mid v)}[h_j]-{\mathrm{\partial }}_{W_{ij}}\log Z$; the same calculation gives ${\mathrm{\partial }}_{a_i}\log p(v)=v_i-{\mathrm{\partial }}_{a_i}\log Z$ and ${\mathrm{\partial }}_{b_j}\log p(v)={\mathbb{E}}_{p(h\mid v)}[h_j]-{\mathrm{\partial }}_{b_j}\log Z$.

The product conditional mass functions in steps 2.1–2.2 make the conditional expectations of any bounded measurable functions of disjoint coordinates factor into the product of their conditional expectations. By [F1]–[F2], each layer is therefore conditionally independent given the other in the stated sigma-algebra sense. This use of the conditional-expectation convention is the only proof point needing [A1].

Differentiating the finite partition sum yields ${\mathrm{\partial }}_{W_{ij}}\log Z={\sum }_{v^{\mathrm{\prime }},h^{\mathrm{\prime }}}{v}_i^{\mathrm{\prime }}{h}_j^{\mathrm{\prime }}p(v^{\mathrm{\prime }},h^{\mathrm{\prime }})={\mathbb{E}}_p[V_i{H}_j]$, ${\mathrm{\partial }}_{a_i}\log Z={\mathbb{E}}_p[V_i]$, and ${\mathrm{\partial }}_{b_j}\log Z={\mathbb{E}}_p[H_j]$. Substitution in step 2.3 gives the exact gradients in the Definition.

For each visible state $v$, ${\sum }_{v^{\mathrm{\prime }}}K(v,v^{\mathrm{\prime }})={\sum }_{h}p(h\mid v){\sum }_{v^{\mathrm{\prime }}}p(v^{\mathrm{\prime }}\mid h)={\sum }_{h}p(h\mid v)=1$. Thus K is a stochastic matrix and each $q_k^v$ is a well-defined finite probability distribution. Replacing the exact negative-phase marginal $p(v^{\mathrm{\prime }})$ in the gradient formulas by $q_k^v(v^{\mathrm{\prime }})$ gives the finite contrastive-divergence statistic defined in the Statement.

Fix distinct $i,j$. By [A1], [F1] and [F2], the joint density of $X$ is the stated nondegenerate Gaussian density; here positive definiteness is an additional hypothesis beyond the singular-covariance multivariate-normal definition.

Reorder coordinates so $Y=(X_i,X_j)$ and $Z=(X_k{)}_{k\mathrm{\notin }\{i,j\}}$. In the correspondingly permuted precision matrix write $\mathrm{\Theta }=\left(\begin{array}{cc}A& B\\ B^T& C\end{array}\right)$, where $A=\left(\begin{array}{cc}a& c\\ c& d\end{array}\right)$. For every nonzero $u\in {\mathbb{R}}^2$, $u^{T}Au=(u,0{)}^T\mathrm{\Theta }(u,0)>0$, so $A$ is positive definite. In particular, $a>0$ and evaluating the quadratic form of $A$ at $(-c\mathrm{/}a,1)$ gives $d-c^2\mathrm{/}a>0$, hence $\det A=ad-c^2>0$. For $p=2$, $Z$ ranges over the singleton ${\mathbb{R}}^0$ and the density of $Y$ is unconditional.

For $p>2$, [F7] supplies the countable-choice hypothesis of [F6], so the joint Euclidean density in [F2] is also a density relative to ${\lambda }_2\otimes {\lambda }_{p-2}$. Put $u=y-{\mu }_Y$, $w=z-{\mu }_Z$, and $D=C-B^T{A}^{-1}B$. For nonzero $w$, $(-A^{-1}Bw,w{)}^T\mathrm{\Theta }(-A^{-1}Bw,w)=w^{T}Dw>0$, hence $D$ is positive definite. Completing the square gives $u^{T}Au+2u^{T}Bw+w^{T}Cw=(u+A^{-1}Bw{)}^{T}A(u+A^{-1}Bw)+w^{T}Dw$. The section integral over $y$ is therefore positive and finite for every $z$: its $y$ factor is the normalized two-dimensional Gaussian density with precision $A$, multiplied by a finite positive constant.

For $p>2$, use the explicit filler probability concentrated at $(0,0)$ on ${\mathbb{R}}^2$ in [F3]; the positive finite marginal sections from step 3.1 show that the conditional law at every $z$ has density $\frac{\sqrt{\det A}}{2\pi }\exp (-\frac{1}{2}(u+A^{-1}Bw{)}^{T}A(u+A^{-1}Bw))$. It is $N_2(m(z),A^{-1})$, where $m(z)={\mu }_Y-A^{-1}B(z-{\mu }_Z)$. For $p=2$, conditioning on the constant $Z$ is unconditional; since $A=\mathrm{\Theta }$ in this case, the law is $N_2(m(\ast ),A^{-1})$ with $m(\ast )={\mu }_Y$. Thus in either case the conditional law is $N_2(m(z),A^{-1})$, interpreting $z=\ast$ when $p=2$.

If ${\mathrm{\Theta }}_{ij}=c=0$, then $A$ is diagonal. For $p>2$, the conditional density in step 4.1 is a product of two univariate Gaussian densities for every $z$; for $p=2$, the unconditional bivariate density is that product. For $p>2$, [F3] gives the conditional-event identities and [F9] identifies the corresponding conditional expectations; for $p=2$, conditional expectations given the constant $Z$ are ordinary expectations. In either case, simple-function approximation extends the indicator identities to bounded measurable functions. Tonelli's theorem [F10], applied to the positive and negative parts of the product integrand, factors the integrals against the product density. Hence the conditional expectations factor for all bounded measurable coordinate functions, and [F4] gives $X_i\perp \text{\hspace{0.17em}⁣}\text{\hspace{0.17em}⁣}\text{\hspace{0.17em}⁣}\perp X_j\mid \sigma (Z)$.

Conversely, suppose $X_i\perp \text{\hspace{0.17em}⁣}\text{\hspace{0.17em}⁣}\text{\hspace{0.17em}⁣}\perp X_j\mid \sigma (Z)$. Apply [F4] to the four bounded pairs formed from $\cos (X_i),\sin (X_i)$ and $\cos (X_j),\sin (X_j)$ and combine the resulting identities. For $p>2$, [F3] and [F9] identify these conditional expectations with integrals under the conditional law in step 4.1, so for almost every $z$ its joint characteristic function at $(1,1)$ equals the product of its conditional marginal characteristic functions. For $p=2$ the same identity is unconditional, at the unique value $z=\ast$. By [F1] and [F5], in either case the joint characteristic function is $\exp (i(m_1(z)+m_2(z))-\frac{1}{2}((A^{-1}{)}_{11}+2(A^{-1}{)}_{12}+(A^{-1}{)}_{22}))$, whereas the product of the marginals omits the middle covariance term. Equality implies $\exp (-(A^{-1}{)}_{12})=1$, and since $(A^{-1}{)}_{12}$ is real this forces $(A^{-1}{)}_{12}=0$.

The inverse of $A$ is $(ad-c^2{)}^{-1}\left(\begin{array}{cc}d& -c\\ -c& a\end{array}\right)$. Since $\det A>0$, $(A^{-1}{)}_{12}=0$ if and only if $c=0$, which is exactly ${\mathrm{\Theta }}_{ij}=0$. Steps 5.1 and 5.2 prove the two implications; by [F8], these are exactly the nonedge and conditional-independence cases of the Gaussian pairwise Markov graph.

Each independent observation with the common law in [F3] has the density in [F2]. On the one-point source $\{\ast \}$, the constant law $K_r(\ast ,\cdot )=N_p(\mu ,\mathrm{\Sigma })$ is a probability kernel: its sections are probability measures and each evaluation map is constant, hence measurable [F4]. The density $f_r(\ast ,x)=f_{\mu ,\mathrm{\Theta }}(x)$ is continuous in $x$ and constant in $\ast$, hence jointly measurable; [F5] gives the product law density ${\prod }_{r=1}^N{f}_{\mu ,\mathrm{\Theta }}(x_r)$. By [F6], independence identifies this product law with the joint observation law. Its log-likelihood for fixed $\mu$ and $\mathrm{\Theta }$ is $L(\mu ,\mathrm{\Theta })=-\frac{Np}{2}\log (2\pi )+\frac{N}{2}\log \det \mathrm{\Theta }-\frac{1}{2}{\sum }_{r=1}^N(x_r-\mu {)}^T\mathrm{\Theta }(x_r-\mu )\mathrm{.}$

Write $y_r=x_r-\bar{x}$, so ${\sum }_r{y}_r=0$ and $x_r-\mu =y_r+(\bar{x}-\mu )$. Expanding and summing removes the cross term: ${\sum }_{r=1}^N(x_r-\mu {)}^T\mathrm{\Theta }(x_r-\mu )={\sum }_{r=1}^N{y}_r^T\mathrm{\Theta }{y}_r+N(\mu -\bar{x}{)}^T\mathrm{\Theta }(\mu -\bar{x})\mathrm{.}$

Positive definiteness of $\mathrm{\Theta }$ makes the last quadratic term nonnegative and equal to zero exactly at $\mu =\bar{x}$. Therefore $\bar{x}$ uniquely maximizes $L(\mu ,\mathrm{\Theta })$ over $\mu$.

By [F1], ${\sum }_r{y}_r{y}_r^T=NS$. Cyclicity of the finite matrix trace gives ${\sum }_r{y}_r^T\mathrm{\Theta }{y}_r={\sum }_r\mathrm{tr}(\mathrm{\Theta }{y}_r{y}_r^T)=N\mathrm{tr}(S\mathrm{\Theta })$. Substituting this identity and $\mu =\bar{x}$ into step 1.1 yields the displayed profile log-likelihood and its normalization.

Initially $W$ is symmetric, and for every nonzero $x\in {\mathbb{R}}^p$, $x^{T}Wx=x^{T}Sx+\lambda \mathrm{\parallel }x{\mathrm{\parallel }}^2\ge \lambda \mathrm{\parallel }x{\mathrm{\parallel }}^2>0$. Thus initialization is positive definite. If the current $W$ is positive definite, reindexing both its rows and columns only reindexes this quadratic form. For a nonzero $z\in {\mathbb{R}}^{p-1}$, extending $z$ by a last zero in the temporary ordering gives $z^{T}Vz>0$. Hence $V$ is positive definite and $V_{\mathrm{\ell }\mathrm{\ell }}=e_{\mathrm{\ell }}^{T}V{e}_{\mathrm{\ell }}>0$ at each coordinate.

Hold all coordinates except ${\beta }_{\mathrm{\ell }}=t$ fixed, and write $v=V_{\mathrm{\ell }\mathrm{\ell }}>0$ and $r=s_{\mathrm{\ell }}-{\sum }_{k\mathrm{\ne }\mathrm{\ell }}{V}_{\mathrm{\ell }k}{\beta }_k$. Expanding the symmetric quadratic gives $Q_{V,s}(\beta )=q(t)+C$, where $q(t)=\frac{1}{2}v{t}^2-rt+\lambda \mathrm{\mid }t\mathrm{\mid }$ and $C$ is independent of $t$. Put $a=(r-\lambda )\mathrm{/}v$ if $r>\lambda$, $a=0$ if $\mathrm{\mid }r\mathrm{\mid }\le \lambda$, and $a=(r+\lambda )\mathrm{/}v$ if $r<-\lambda$. In these three cases respectively put $u=1$, $u=r\mathrm{/}\lambda$, and $u=-1$. Then $va-r+\lambda u=0$, $\mathrm{\mid }u\mathrm{\mid }\le 1$, and $ua=\mathrm{\mid }a\mathrm{\mid }$, so $\mathrm{\mid }t\mathrm{\mid }-\mathrm{\mid }a\mathrm{\mid }-u(t-a)=\mathrm{\mid }t\mathrm{\mid }-ut\ge 0$. Consequently $q(t)-q(a)=\frac{1}{2}v(t-a{)}^2+\lambda (\mathrm{\mid }t\mathrm{\mid }-\mathrm{\mid }a\mathrm{\mid }-u(t-a))\ge \frac{1}{2}v(t-a{)}^2$, and the lower bound is strictly positive whenever $t\mathrm{\ne }a$. Thus $a=\mathsf{S}\mathsf{o}\mathsf{f}\mathsf{t}(r,\lambda )\mathrm{/}v$ is the unique coordinate minimizer. All quantities are real and division is by a positive number, so finitely many such updates produce a finite real vector $\beta$; the nonpositive-diagonal guard cannot fire while the current $W$ is positive definite.

For such a block, $V\beta$ and $d$ are finite and the block candidate $\left(\begin{array}{cc}V& V\beta \\ (V\beta {)}^T& d\end{array}\right)$ is symmetric. Restoring original indexing preserves symmetry. By [F2] any symmetric candidate $A$ has a representation $A=UD{U}^T$ with $U$ orthogonal and $D$ real diagonal. The identity $x^{T}Ax={\sum }_i{D}_{ii}(U^{T}x{)}_i^2$ shows that $A$ is positive definite if all $D_{ii}>0$; conversely applying positive definiteness to $x=U{e}_i$ shows that every $D_{ii}>0$. Thus the candidate guard accepts exactly the positive-definite candidates. Initialization, the within-block reasoning and this acceptance rule give induction over all reached blocks: every accepted $W$ is symmetric positive definite, and each later block again has the properties proved above until a candidate guard returns failure.

There are at most $Rp$ blocks, each with $K(p-1)$ coordinate updates and one finite-dimensional eigenvalue test; early failure only decreases these counts. Every reached operation is well-defined by the preceding steps and the stipulated exact matrix operations, so execution terminates with failure or completes all blocks. In the latter case the final $W=UD{U}^T$ has positive diagonal eigenvalues by step 3.1. Hence $W^{-1}=U{D}^{-1}{U}^T$ exists, is symmetric and has $x^T{W}^{-1}x={\sum }_i{D}_{ii}^{-1}(U^{T}x{)}_i^2>0$ for every nonzero $x$. This proves every asserted property without a convergence or global-optimizer inference.

By the all-entry objective definition [F1] and the nonnegative-operator factorization [F2], write $S=R^{T}R$. Thus $\mathrm{tr}(S\mathrm{\Theta })=\mathrm{tr}(R\mathrm{\Theta }{R}^T)\ge 0$ for $\mathrm{\Theta }\succ 0$. Also $\mathrm{\parallel }\mathrm{\Theta }{\mathrm{\parallel }}_1={\sum }_{i,j}\mathrm{\mid }{\mathrm{\Theta }}_{ij}\mathrm{\mid }\ge {\sum }_i{\mathrm{\Theta }}_{ii}=\mathrm{tr}\mathrm{\Theta }$, since positive-definite matrices have positive diagonal entries.

The cone ${\mathbb{S}}_{++}^p$ is open and convex: positive quadratic forms remain positive under sufficiently small symmetric perturbations, and convex combinations of positive quadratic forms remain positive. For symmetric $H$, determinant expansion gives ${\frac{d}{dt}\log \det (\mathrm{\Theta }+tH)\mid }_{t=0}=\mathrm{tr}({\mathrm{\Theta }}^{-1}H)$. Differentiating $(\mathrm{\Theta }+tH)(\mathrm{\Theta }+tH{)}^{-1}=I$ gives ${\frac{d}{dt}(\mathrm{\Theta }+tH{)}^{-1}\mid }_{t=0}=-{\mathrm{\Theta }}^{-1}H{\mathrm{\Theta }}^{-1}$, so the second directional derivative is $-\mathrm{tr}({\mathrm{\Theta }}^{-1}H{\mathrm{\Theta }}^{-1}H)$.

By [F3], the eigenvalues $t_1,\dots ,t_p$ of $\mathrm{\Theta }\succ 0$ are positive and $\log \det \mathrm{\Theta }={\sum }_k\log t_k$, $\mathrm{tr}\mathrm{\Theta }={\sum }_k{t}_k$. Hence $F_{S,\lambda }(\mathrm{\Theta })\le {\sum }_{k=1}^{p}h(t_k)$ for $h(t)=\log t-\lambda t$. The function $h$ is bounded above and tends to $-\mathrm{\infty }$ as $t\downarrow 0$ or $t\to \mathrm{\infty }$. For any $c$, the upper level set $F_{S,\lambda }(\mathrm{\Theta })\ge c$ therefore forces each eigenvalue into a compact interval $[a_c,b_c]\subset (0,\mathrm{\infty })$: if $M={\sup }_{t>0}h(t)$, then each $h(t_k)\ge c-(p-1)M$. These matrices form a closed bounded set of symmetric matrices contained in ${\mathbb{S}}_{++}^p$, hence a compact set. Taking $c=F_{S,\lambda }(I)$ gives a nonempty compact upper level set, and continuity of $F$ yields a maximizer.

Let $B={\mathrm{\Theta }}^{-1\mathrm{/}2}H{\mathrm{\Theta }}^{-1\mathrm{/}2}$, using the positive square root from the orthonormal eigenbasis [F3]. Then $B$ is symmetric and $\mathrm{tr}({\mathrm{\Theta }}^{-1}H{\mathrm{\Theta }}^{-1}H)=\mathrm{tr}(B^2)=\mathrm{\parallel }B{\mathrm{\parallel }}_F^2$, strictly positive when $H\mathrm{\ne }0$. Thus the second derivative of $\log \det$ along every nonconstant segment is strictly negative, giving its strict chord inequality. The Hessian of $-\log \det$ is positive semidefinite by step 1.2, so [F4] gives convexity in the chord sense of [F5]. The trace term is linear and the entrywise ${\mathrm{\ell }}^1$ norm is convex by the triangle inequality, so its negative is concave; adding these terms preserves the strict chord inequality. Hence $F_{S,\lambda }$ is strictly concave on its open convex domain.

Put $G={\mathrm{\Theta }}^{-1}-S$. At a maximizer every feasible one-coordinate directional derivative is at most zero, and when the objective is differentiable in both directions its derivative is zero. Since ${\mathrm{\Theta }}_{ii}>0$, the diagonal penalty is locally differentiable with slope $1$, so the diagonal equation is $G_{ii}-\lambda =0$ and ${\mathrm{\Gamma }}_{ii}=1$. For $i<j$, a symmetric perturbation by $t(E_{ij}+E_{ji})$ has smooth derivative $2G_{ij}t$ and changes the penalty by $2\lambda (\mathrm{\mid }{\mathrm{\Theta }}_{ij}+t\mathrm{\mid }-\mathrm{\mid }{\mathrm{\Theta }}_{ij}\mathrm{\mid })$. If ${\mathrm{\Theta }}_{ij}\mathrm{\ne }0$, testing both signs of small $t$ gives $G_{ij}=\lambda \mathrm{sign}({\mathrm{\Theta }}_{ij})$. If ${\mathrm{\Theta }}_{ij}=0$, taking $t>0$ gives $G_{ij}\le \lambda$ and taking $t<0$ gives $G_{ij}\ge -\lambda$. Defining ${\mathrm{\Gamma }}_{ij}={\mathrm{\Gamma }}_{ji}=G_{ij}\mathrm{/}\lambda$ at zero and the sign otherwise yields the stated KKT equation for the all-entry penalty [F1].

Conversely, suppose $\mathrm{\Theta }\succ 0$ and a symmetric $\mathrm{\Gamma }$ satisfies the stated entrywise conditions and ${\mathrm{\Theta }}^{-1}-S-\lambda \mathrm{\Gamma }=0$. Put $T={\mathrm{\Theta }}^{-1\mathrm{/}2}{\mathrm{\Theta }}^{\mathrm{\prime }}{\mathrm{\Theta }}^{-1\mathrm{/}2}$ for any ${\mathrm{\Theta }}^{\mathrm{\prime }}\succ 0$; its eigenvalues $t_r$ are positive. Since $\log t\le t-1$ for $t>0$ (the function $t-1-\log t$ has its minimum $0$ at $t=1$), $\log \det {\mathrm{\Theta }}^{\mathrm{\prime }}-\log \det \mathrm{\Theta }={\sum }_r\log t_r\le \mathrm{tr}(T)-p=\mathrm{tr}({\mathrm{\Theta }}^{-1}({\mathrm{\Theta }}^{\mathrm{\prime }}-\mathrm{\Theta }))$. For each scalar $x$ and each allowed subgradient $\gamma$ of $\mathrm{\mid }x\mathrm{\mid }$, $\mathrm{\mid }y\mathrm{\mid }\ge \mathrm{\mid }x\mathrm{\mid }+\gamma (y-x)$: if $x\mathrm{\ne }0$, this follows from $\mathrm{\mid }y\mathrm{\mid }\ge \mathrm{sign}(x)y$, while for $x=0$ it follows from $\mathrm{\mid }\gamma \mathrm{\mid }\le 1$. Summing gives $\mathrm{\parallel }{\mathrm{\Theta }}^{\mathrm{\prime }}{\mathrm{\parallel }}_1\ge \mathrm{\parallel }\mathrm{\Theta }{\mathrm{\parallel }}_1+{\sum }_{i,j}{\mathrm{\Gamma }}_{ij}({\mathrm{\Theta }}_{ij}^{\mathrm{\prime }}-{\mathrm{\Theta }}_{ij})$ for the all-entry objective [F1]. Combining these bounds with the linear trace term and the KKT equation yields $F_{S,\lambda }({\mathrm{\Theta }}^{\mathrm{\prime }})-F_{S,\lambda }(\mathrm{\Theta })\le \mathrm{tr}(({\mathrm{\Theta }}^{-1}-S-\lambda \mathrm{\Gamma })({\mathrm{\Theta }}^{\mathrm{\prime }}-\mathrm{\Theta }))=0$. Thus every such KKT point is a maximizer. The maximizer exists by step 2.1 and is unique by strict concavity from step 2.2, so the KKT condition is both necessary and sufficient for that unique optimum.