Structured regularization for conditional Gaussian graphical models

Abstract

Conditional Gaussian graphical models are a reparametrization of the multivariate linear regression model which explicitly exhibits (i) the partial covariances between the predictors and the responses, and (ii) the partial covariances between the responses themselves. Such models are particularly suitable for interpretability since partial covariances describe direct relationships between variables. In this framework, we propose a regularization scheme to enhance the learning strategy of the model by driving the selection of the relevant input features by prior structural information. It comes with an efficient alternating optimization procedure which is guaranteed to converge to the global minimum. On top of showing competitive performance on artificial and real datasets, our method demonstrates capabilities for fine interpretation, as illustrated on three high-dimensional datasets from spectroscopy, genetics, and genomics.

Appendix: Proofs

Most of the proofs rely on basic algebra and properties of the trace and the $\mathrm {vec}$ operators (see e.g. Harville 1997), in particular

$$\begin{aligned}&\mathrm {tr}(\mathbf {A}^T \mathbf {B}\mathbf {C}\mathbf {D}^T) = \mathrm {vec}(\mathbf {A})^T (\mathbf {D}\otimes \mathbf {B}) \mathrm {vec}(\mathbf {C}),\\&\mathrm {vec}(\mathbf {A}\mathbf {B}\mathbf {C}) = (\mathbf {C}^T \otimes \mathbf {A}) \mathrm {vec}(\mathbf {B}), \\&(\mathbf {A}\otimes \mathbf {B}) (\mathbf {C}\otimes \mathbf {D}) = (\mathbf {A}\mathbf {C}) \otimes (\mathbf {C}\mathbf {D}). \end{aligned}$$

1.1 Derivation of Proposition 1

Concerning the two regularization terms in (7), we have $\Vert {\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}}\Vert _1 = \Vert {\varvec{\omega }}\Vert _1$ since the $\ell _1$-norm applies element-wise here and

$$\begin{aligned} \mathrm {tr}({\varvec{\varOmega }}_{\mathbf {y}\mathbf {x}} \mathbf {L}{\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}}{\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}}^{-1}) = {\varvec{\omega }}^T ({\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}}^{-1} \otimes \mathbf {L}) {\varvec{\omega }}. \end{aligned}$$

As for the log-likelihood (5), we work on the trace term:

$$\begin{aligned}&\mathrm {tr}\left( \left( \mathbf {Y}+ \mathbf {X}{\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}}{\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}}^{-1} \right) ^T \left( \mathbf {Y}+ \mathbf {X}{\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}}{\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}}^{-1} \right) {\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}}\right) \\&\quad = \mathrm {vec}\left( \mathbf {Y}+ \mathbf {X}{\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}}{\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}}^{-1}\right) ^T \left( {\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}}\otimes \mathbf {I}_n\right) \mathrm {vec}\left( \mathbf {Y}+ \mathbf {X}{\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}}{\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}}^{-1}\right) \\&\quad = \left\| \left( {\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}}\otimes \mathbf {I}_n\right) ^{1/2} \left( \mathrm {vec}(\mathbf {Y}) + \left( {\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}}^{-1} \otimes \mathbf {X}\right) \mathrm {vec}({\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}}) \right) \right\| _2^2 \\&\quad = \left\| \mathrm {vec}(\mathbf {Y}{\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}}^{1/2}) + \left( {\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}}^{-1/2}\otimes \mathbf {X}\right) {\varvec{\omega }}\right\| _2^2. \end{aligned}$$

The rest of the proof is straightforward.

1.2 Convexity lemma

Lemma 1

The function

$$\begin{aligned} -\frac{1}{n}\log L({\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}},{\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}}) + \frac{\lambda _2}{2} \mathrm {tr}\left( {\varvec{\varOmega }}_{\mathbf {y}\mathbf {x}}\mathbf {L}{\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}}{\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}}^{-1}\right) \end{aligned}$$

is jointly convex in $({\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}},{\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}})$ and admits at least one global minimum which is unique when $n\ge q$ and $(\lambda _2\mathbf {L}+ \mathbf {S}_{\mathbf {x}\mathbf {x}})$ is positive definite.

The convexity of $-\frac{1}{n}\log L({\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}},{\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}})$ is proved in Yuan and Zhang (2014) (Proposition 1). Similar arguments can be straightforwardly applied in the case at hand. Existence of the global minimum is related to strict convexity in both ${\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}}$ and ${\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}}$, where direct differentiation leads to the corresponding conditions.

1.3 Proof of Proposition 2

The convexity of criterion $J({\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}},{\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}})$ in $({\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}},{\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}})$ is straightforward thanks to Lemma 1 and considering the fact that $\Vert {\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}}\Vert _1$ is also convex. One can then apply the results developed in Tseng and Yun (2009), Tseng (2001) on the convergence of block coordinate descent for the minimization of nonsmooth separable function. Since (7) is clearly separable in $({\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}}, {\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}})$ for the nonsmooth part induced by the $\ell _1$-norm, the alternating scheme is guaranteed to converge to the unique global minimum under the assumption of Lemma 1. It remains to show that the two convex optimization subproblems (9a) and (9b) can be (efficiently) solved in practice.

Firstly, (9b) can be recast as an Elastic-Net problem, which in turn can be recast as a LASSO problem (see, e.g. Zou and Hastie 2005; Slawski et al. 2010). This is straightforward thanks to Proposition 1: when $\hat{{\varvec{\varOmega }}}_{\mathbf {y}\mathbf {y}}$ is fixed, solution to (9b) can be obtained via

(13)

where $\mathbf {A},\mathbf {b}$ and $\tilde{\mathbf {L}}$ are defined by

$$\begin{aligned} \mathbf {A}= & {} \left( \hat{{\varvec{\varOmega }}}_{\mathbf {y}\mathbf {y}}^{-1/2} \otimes \mathbf {X}/\sqrt{n} \right) ,\\ \mathbf {b}= & {} - \mathrm {vec}\left( \mathbf {Y}\hat{{\varvec{\varOmega }}}_{\mathbf {y}\mathbf {y}}^{1/2}\right) / \sqrt{n} \text { and } \tilde{\mathbf {L}} = \hat{{\varvec{\varOmega }}}_{\mathbf {y}\mathbf {y}}^{-1} \otimes \mathbf {L}. \end{aligned}$$

Secondly, we can solve analytically (9a) with simple matrix algebra. By differentiation of the objective (7) over ${\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}}^{-1}$ we obtain the quadratic form

$$\begin{aligned} {\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}} \mathbf {S}_{\mathbf {y}\mathbf {y}} {\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}} - {\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}} = {\varvec{\varOmega }}_{\mathbf {y}\mathbf {x}}(\lambda _2 \mathbf {L}+ \mathbf {S}_{\mathbf {x}\mathbf {x}}){\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}}. \end{aligned}$$

After right multiplying both sides by $\mathbf {S}_{\mathbf {y}\mathbf {y}}$, it becomes obvious that ${\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}} \mathbf {S}_{\mathbf {y}\mathbf {y}}$ and ${\varvec{\varOmega }}_{\mathbf {y}\mathbf {x}}(\lambda _2 \mathbf {L}+ \mathbf {S}_{\mathbf {x}\mathbf {x}}){\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}}\mathbf {S}_{\mathbf {y}\mathbf {y}}$ commute and thus share the same eigenvectors $\mathbf {U}$. Besides, it induces the relationship $\eta _j^2 -\eta _j = \zeta _j$ between their respective eigenvalues $\eta _j$ and $\zeta _j$, and we are looking for the positive solution of $\eta _j$. To do so, first note that we may assume that ${\varvec{\varOmega }}_{\mathbf {y}\mathbf {x}} (\lambda _2 \mathbf {L}+ \mathbf {S}_{\mathbf {x}\mathbf {x}}) {\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}}$ and $\mathbf {S}_{\mathbf {y}\mathbf {y}}$ are positive definite, when ${\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}} \ne \mathbf {0}$ and $n\ge q$; and second, recall that if a matrix is the product of two positive definite matrices then its eigenvalues are positive. Hence, $\zeta _j>0$ and the positive solution of $\eta _j$ is $\eta _j = (1 + \sqrt{1+4\zeta _j})/2$. We thus obtain

$$\begin{aligned} \hat{{\varvec{\varOmega }}}_{\mathbf {y}\mathbf {y}} = \mathbf {U}\mathrm {diag}(\varvec{\eta }) \mathbf {U}^{-1} \mathbf {S}_{\mathbf {y}\mathbf {y}}^{-1}. \end{aligned}$$

(14)

Direct inversion yields

$$\begin{aligned} \hat{{\varvec{\varOmega }}}_{\mathbf {y}\mathbf {y}} = \mathbf {U}\mathrm {diag}(\varvec{\eta }/\varvec{\zeta }) \mathbf {U}^{-1} \hat{{\varvec{\varOmega }}}_{\mathbf {y}\mathbf {x}} (\lambda _2\mathbf {L}+ \mathbf {S}_{\mathbf {x}\mathbf {x}}) \hat{{\varvec{\varOmega }}}_{\mathbf {x}\mathbf {y}} (=\hat{\mathbf {R}}^{-1}), \end{aligned}$$

To get an expression for $\hat{{\varvec{\varOmega }}}_{\mathbf {y}\mathbf {y}}$ which does not require additional matrix inversion, just note that

$$\begin{aligned} \mathbf {S}_{\mathbf {y}\mathbf {y}}^{-1}= & {} \left( {\varvec{\varOmega }}_{\mathbf {y}\mathbf {x}}\hat{{\varvec{\varSigma }}}_{\mathbf {x}\mathbf {x}}^{\lambda _2}{\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}} \mathbf {S}_{\mathbf {y}\mathbf {y}}\right) ^{-1} {\varvec{\varOmega }}_{\mathbf {y}\mathbf {x}}\hat{{\varvec{\varSigma }}}_{\mathbf {x}\mathbf {x}}^{\lambda _2}{\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}} \\= & {} \mathbf {U}\mathrm {diag}(\varvec{\zeta }^{-1})\mathbf {U}^{-1}{\varvec{\varOmega }}_{\mathbf {y}\mathbf {x}} \hat{{\varvec{\varSigma }}}_{\mathbf {x}\mathbf {x}}^{\lambda _2}{\varvec{\varOmega }}_{\mathbf {x}\mathbf {y}} \end{aligned}$$

where $\hat{{\varvec{\varSigma }}}_{\mathbf {x}\mathbf {x}}^{\lambda _2} = (\lambda _2\mathbf {L}+ \mathbf {S}_{\mathbf {x}\mathbf {x}})$. Combined with (14), this last equality leads to

$$\begin{aligned} \hat{{\varvec{\varOmega }}}_{\mathbf {y}\mathbf {y}}^{-1} = \mathbf {S}_{\mathbf {y}\mathbf {y}}\mathbf {U}\mathrm {diag}(\varvec{\eta }^{-1}) \mathbf {U}^{-1} (=\hat{\mathbf {R}}). \end{aligned}$$

Finally, in the particular case where $\hat{{\varvec{\varOmega }}}_{\mathbf {x}\mathbf {y}} = \mathbf {0}$, $\hat{{\varvec{\varOmega }}}_{\mathbf {y}\mathbf {y}} = \mathbf {S}_{\mathbf {y}\mathbf {y}}^{-1}$.

1.4 Derivation of Proposition 3

To apply the results developed in Tibshirani and Taylor (2012) that rely on the well-known Stein’s Lemma (Stein 1981), we basically need to recast our problem as a classical LASSO applied on a Gaussian linear regression model such that the response vector follows a normal distribution of the form $\mathscr {N}({\varvec{\mu }},\sigma \mathbf {I})$. This is straightforward by means of Proposition 1: in the same way as we derive expression (13), we can go a little further and reach the following LASSO formulation

(15)

with $\mathbf {A},\mathbf {b}$ and $\tilde{\mathbf {L}}$ defined as in (13). From model (4), it is not difficult to see that $\mathbf {b}$ corresponds to an uncorrelated vector form of $\mathbf {Y}$ so as

$$\begin{aligned} \begin{pmatrix} \mathbf {b} \\ \mathbf {0}\end{pmatrix} = \begin{pmatrix} -\mathrm {vec}(\mathbf {Y}{\varvec{\varOmega }}_{\mathbf {y}\mathbf {y}}^{1/2}) \\ \mathbf {0}\end{pmatrix} \sim \mathscr {N}({\varvec{\mu }}, \mathbf {I}_{nq}). \end{aligned}$$

The explicit form of ${\varvec{\mu }}$ is of no interest here. The point is essentially to underline that the response vector is uncorrelated in (15), which allows us to apply Theorem 2 of Tibshirani and Taylor (2012) and results therein, notably for the Elastic-Net. By these means, an unbiased estimator of the degrees of freedom in (15) can be written as a function of the active set $\mathscr {A}$ in $\hat{{\varvec{\omega }}}^{\lambda _1,\lambda _2}$:

$$\begin{aligned} \hat{\mathrm {df}}_{\lambda _1,\lambda _2} = \mathrm {tr}\left( \mathbf {A}_{\mathscr {A}} \left( (\mathbf {A}^T\mathbf {A}+ \lambda _2 \tilde{\mathbf {L}})_{\mathscr {A}\mathscr {A}}\right) ^{-1} \mathbf {A}_{\mathscr {A}}^T \right) . \end{aligned}$$

Routine simplifications lead to the desired result for the degrees of freedom.