Empirical Bayesian Learning in AR Graphical Models

TL;DR

This paper introduces an empirical Bayes approach for learning sparse autoregressive graphical models, demonstrating its effectiveness through numerical experiments in high-dimensional settings.

Contribution

It proposes a novel empirical Bayes estimator for sparse and latent-variable autoregressive graphical models, advancing Bayesian methods in this domain.

Findings

Effective in high-dimensional settings

Improves learning of sparse graphical models

Numerical experiments confirm benefits of Bayesian approach

Abstract

We address the problem of learning graphical models which correspond to high dimensional autoregressive stationary stochastic processes. A graphical model describes the conditional dependence relations among the components of a stochastic process and represents an important tool in many fields. We propose an empirical Bayes estimator of sparse autoregressive graphical models and latent-variable autoregressive graphical models. Numerical experiments show the benefit to take this Bayesian perspective for learning these types of graphical models.

Figures

Equations

$${\mathcal{Q}}_{m,n}=\{\Delta X{\Delta }^{\ast }s.t.X\in {\mathbf{Q}}_{m(n+1)}\}$$

$$y(t)=-\sum _{k=1}^n{A}_{k}y(t-k)+e(t).$$

$$\Phi (e^{i\vartheta })=\sum _{k\in \mathbb{Z}}{e}^{-ik\vartheta }{R}_k,\vartheta \in [-\pi ,\pi ]$$

$$(j,h)\notin \mathcal{E}⟺{\chi }_{\{j\}}\perp {\chi }_{\{h\}}|{\chi }_{\mathcal{V}\setminus \{j,h\}},j\ne h.$$

$$\hat{\Sigma }=$$

$$s.t.\; \Sigma \succ 0.$$

$$\ell (y^N;\Sigma )=-\frac{N-n}{2}\int \log |\Sigma |+\frac{N-n}{2}⟨{\Phi }_{y^N},\Sigma ⟩+c$$

$${\Phi }_{y^N}={\hat{R}}_0+\frac{1}{2}\sum _{k=1}^n{\hat{R}}_k{e}^{-ik\vartheta }+{\hat{R}}_k^T{e}^{ik\vartheta }$$

$$\hat{R}=[{\hat{R}}_0\dots {\hat{R}}_n],{\hat{R}}_k=\frac{1}{N-n}\sum _{t=1}^{N-k}y(t+k)y(t{)}^T$$

$${\bar{h}}_{\infty }$$

$$q_{jh}(\Sigma ):=\max \{|(S_0{)}_{jh}|,\underset{k=1\dots n}{\max }|(S_k{)}_{jh}|,\underset{k=1\dots n}{\max }|(S_k{)}_{hj}|\}$$

$$\hat{X}=$$

$$s.t.X_{00}\succ 0,X⪰0$$

$$\hat{\Sigma }=$$

$$s.t.\; \Sigma \succ 0.$$

$$\hat{X}=$$

$$s.t.X_{00}\succ 0,X⪰0$$

$$\underset{\begin{array}{c}Z\in {\mathbf{M}}_{m,n}\\ W\in {\mathbf{Q}}_m\end{array}}{\max }\log |W|+m$$

$$s.t.W\succ 0$$

$$\displaystyle\hskip 19.91684pt\mathrm{T}(\hat{R})+\mathrm{T}(Z)\succeq\left[\begin{array}[]{cc}W&0\\ 0&0\end{array}\right]$$

$$\sum _{k=0}^n|(Z_k{)}_{jh}|+|(Z_k{)}_{hj}|\le \frac{2{\gamma }_{jh}}{N-n},j>h$$

$$\sum _{k=0}^n|(Z_k{)}_{jj}|\le \frac{2{\gamma }_{jj}}{N-n},j=h.$$

$$\underset{\begin{array}{c}X\in {\mathbf{Q}}_{m(n+1)}\end{array}}{\mathrm{argmin}}\sum _{j\ge h}\log (q_{jh}(\mathrm{D}(X))+{\epsilon }_S)+\ell (y^N;X)$$

$$s.t.X_{00}\succ 0,X⪰0.$$

$$\underset{\begin{array}{c}X\in {\mathbf{Q}}_{m(n+1)}\\ u_{jh}\in \mathbb{R},j\ge h\end{array}}{\mathrm{argmin}}\sum _{j\ge h}\log (u_{jh}+{\epsilon }_S)+\ell (y^N;X)$$

$$s.t.X_{00}\succ 0,X⪰0$$

$$u_{jh}\ge q_{jh}(\mathrm{D}(X)).$$

$$f(u,X)=\sum _{j\ge h}\log (u_{jh}+{\epsilon }_S)+\ell (y^N;X).$$

$$g(u$$

$$:=f({\hat{u}}^{(l)},{\hat{X}}^{(l)})+⟨{\nabla }_{u}f({\hat{u}}^{(l)}),u-{\hat{u}}^{(l)}⟩+\ell (y^N;X)$$

$$=f({\hat{u}}^{(l)},{\hat{X}}^{(l)})+\sum _{j\ge h}\frac{u_{jh}-{\hat{u}}_{jh}^{(l)}}{{\hat{u}}_{jh}^{(l)}+{\epsilon }_S}+\ell (y^N;X).$$

$$({\hat{u}}^{(l+1)},{\hat{X}}^{(l+1)})=$$

$$s.t.X_{00}\succ 0,X⪰0$$

$$u_{jh}\ge q_{jh}(\mathrm{D}(X)).$$

$$(u^{(l+1)},X^{(l+1)})=$$

$$s.t.X_{00}\succ 0,X⪰0$$

$$u_{jh}\ge q_{jh}(\mathrm{D}(X)).$$

$${\hat{X}}^{(l+1)}=$$

$$s.t.X_{00}\succ 0,X⪰0.$$

$${\hat{\gamma }}_{jh}^{(l+1)}:=\frac{1}{q_{jh}(\mathrm{D}({\hat{X}}^{(l)}))+{\epsilon }_S}$$

Full text

1 Introduction

In modern applications many variables are accessible to observation. In some cases the latter can be modeled with a high dimensional Gaussian random vector. To gain some insight about the relation among those variables we can attach to it a graphical model (Lauritzen, 1996; Willsky, 2002). The latter is an undirected graph wherein nodes correspond to the components (i.e. variables) of the random vector and there is the lack of an edge between two nodes if and only if the corresponding variables are conditionally independent given the others. It turns out that sparse graphical models, i.e. graphs with few edges, have the inverse covariance matrix of the random vector which is sparse. Then, the problem of estimating a sparse graphical model from the observed data can be formulated as a regularization problem: find such an inverse matrix which minimizes the negative log-likelihood and a regularization term inducing sparsity (Banerjee et al., 2008).

An important aspect is that variables are typically measured over time and can thus be modeled as a high dimensional autoregressive (AR) Gaussian stationary stochastic process. Then, we can attach a graphical model describing the conditional dependence relations among the variables. It is possible to prove that sparse graphical models have the inverse power spectral density (PSD) of the process which is sparse. Songsiri & Vandenberghe (2010) proposed a regularized estimator for estimating sparse AR graphical models in the same spirit of (Banerjee et al., 2008). Avventi et al. (2013) showed that the aforementioned estimator is a relaxed version of a maximum entropy estimator. The latter solves a covariance extension problem whose dual problem does coincide with the one proposed in Songsiri et al. (2010). As a consequence, the estimator proposed by Songsiri & Vandenberghe (2010) is strictly connected with the generalized moment problems in the sense of Byrnes-Georgiou-Lindquist which have been extensively studied by many researchers, e.g. Byrnes et al. (2000); Ferrante et al. (2012); Karlsson & Georgiou (2013); Zorzi (2015, 2014a, 2015, 2014b). Since then, many other extensions has been proposed: Maanan et al. (2017) proposed a two stage approach to estimate sparse AR graphical models; Alpago et al. (2018) proposed a regularized estimator for sparse graphical models of reciprocal processes; Chandrasekaran et al. (2010), Zorzi & Sepulchre (2016), Maanan et al. (2018), Liégeois et al. (2015) and Ciccone et al. (2018) proposed regularized estimators for the so called latent-variable graphical models.

The regularizers for inducing sparsity in such graphical models are $\ell_{1}$-like norms. The $\ell_{1}$ norm, however, penalizes differently the nonnull coefficients: larger coefficients are penalized more heavily than smaller coefficients. This imbalance produces an estimator with a remarkable mean squared error. Candès et al. (2008) proposed to reduce this imbalance by considering a weighted $\ell_{1}$ norm wherein the weights are computed in an iterative fashion. The resulting procedure is similar to the adaptive lasso and is typically called iterative reweighted algorithm (Wipf & Nagarajan, 2010).

The present paper proposes a regularized estimator in the spirit of Songsiri & Vandenberghe (2010) for sparse AR graphical models where the $\ell_{1}$-like norm is substituted by a weighted $\ell_{1}$-like norm leading to an iterative reweighted procedure. Interestingly, drawing inspiration by Asadi et al. (2009); Scheinberg et al. (2010), the proposed method can be understood as an empirical Bayes approach which provides a suitable updating rule for the weights. Such idea is then extended to the regularized estimator for latent-variable AR graphical models proposed in Zorzi & Sepulchre (2016).

The outline of the paper is as follows. In Section 2 we introduce the problem of estimating sparse AR graphical models. In Section 3 we propose an iterative reweighted for solving such a problem, while in Section 4 we derive the estimator using a Bayesian perspective. Section 5 regards the identification of latent-variable AR graphical models. Section 6 contains some numerical experiments to test the performance of the proposed estimators. Finally, the conclusions are drawn in Section 7.

Notation

The vector space $\mathbb{R}^{m}$ is endowed with the inner product $\left<x,y\right>=x^{T}y$. $x\geq 0$ ($x>0$) with $x\in\mathbb{R}^{m}$ means that all the entries of the vector are nonnegative (positive). The vector space $\mathbb{R}^{m\times m}$ is endowed with the inner product $\left<X,L\right>=\mathop{\rm tr}(XL^{T})$. $\mathbf{Q}_{m}$ denotes the vector space of symmetric matrices of dimension $m\times m$, if $X\in\mathbf{Q}_{m}$ is positive definite (semi-definite) we write $X\succ 0$ ($X\succeq 0$). $|X|$ denotes the determinant of matrix $X\in\mathbf{Q}_{m}$. $\mathrm{diag}(X)$ is the diagonal matrix whose main diagonal coincides with the one of $X$. $(X)_{jh}$ denotes the entry in position $(j,h)$ of matrix $X$. A matrix $A\in\mathbb{R}^{l\times m(n+1)}$ with $l\leq m$ will be partitioned as $A=[\,A_{0}\,A_{1}\,\ldots\,A_{n}\,]$ with $A_{j}\in\mathbb{R}^{l\times m}$. $\mathbf{M}_{m,n}$ is the vector space of matrices $Y:=[\,Y_{0}\,Y_{1}\,\ldots\,Y_{n}\,]$ with $Y_{0}\in\mathbf{Q}_{m}$ and $Y_{1}\ldots Y_{n}\in\mathbb{R}^{m\times m}$. The corresponding inner product is $\left<Y,Z\right>=\mathop{\rm tr}(YZ^{T})$. The linear operator $\mathrm{T}:\mathbf{M}_{m,n}\rightarrow\mathbf{Q}_{m(n+1)}$ constructs a symmetric Toeplitz matrix $\mathrm{T}(Y)$ where its first block row is $[\,Y_{0}\,Y_{1}\,\ldots\,Y_{n}\,].$ The adjoint operator of $\mathrm{T}$ is denoted by $\mathrm{D}:\mathbf{Q}_{m(n+1)}\rightarrow\mathbf{M}_{m,n}$ and defined as follows. If $X\in\mathbf{Q}_{m(n+1)}$ is partitioned as a $n+1\times n+1$ block matrix with $X_{hj}$, $j,h=0\ldots n$, the block in position $(h,j)$ then $\mathrm{D}(X)=[\,\mathrm{D}_{0}(X)\,\ldots\,\mathrm{D}_{n}(X)\,]$ where $\mathrm{D}_{0}(X)=\sum_{h=0}^{n}X_{hh}$, $\mathrm{D}_{k}(X)=2\sum_{h=0}^{n-k}X_{h\;h+k}$, $k=1\ldots n$. We define the index set $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ with $\mathcal{V}:=\{1,2,\ldots m\}$. Functions on the unit circle $\{e^{i\vartheta}\hbox{ s.t. }\vartheta\in[-\pi,\pi]\}$ will be denoted by capital Greek letters, e.g. $\Phi(e^{i\vartheta})$ with $\vartheta\in[-\pi,\pi]$, and the dependence upon $\vartheta$ will be dropped if not needed, e.g. $\Phi$ instead of $\Phi(e^{i\vartheta})$. $\mathrm{L}_{2}^{m\times m}$ denotes the space of $\mathbb{C}^{m\times m}$-valued functions defined on the unit circle which are square integrable. Given $\Phi\in\mathrm{L}^{m\times m}_{2}$, the shorthand notation $\int\Phi$ denotes the integration of $\Phi$ taking place on the unit circle with respect to the normalized Lebesgue measure. Then, the inner product in $\mathrm{L}_{2}^{m\times m}$ is $\left<\Phi,\Sigma\right>=\mathop{\rm tr}\int\Phi\Sigma^{*}$. Given an analytic function $\Lambda\in\mathrm{L}_{2}^{m\times m}$, its (normal) rank is denoted by $\mathrm{rank}(\Lambda)$. If $\Phi(e^{i\vartheta})$ is positive definite (semi-definite) for each $\vartheta\in[-\pi,\pi]$, we will write $\Phi\succ 0$ ($\Phi\succeq 0$). We define the following family of matrix pseudo-polynomials

[TABLE]

where $\Delta(e^{i\vartheta}):=[\,I_{m}\,e^{i\vartheta}I_{m}\,\ldots\,e^{in\vartheta}I_{m}\,]$ is the shift operator. Given a $m$-dimensional stochastic process $y=\{\,y(t),\;t\in\mathbb{Z}\,\}$, $y_{j}$ denotes the $j$-th component (i.e. variable) of $y$. With some abuse of notation, $y(t)$ will both denote a random vector and its sample value. Given a function $f(u,v)$, $\nabla_{u}f(\bar{u},\bar{v})$ and $\nabla_{v}f(\bar{u},\bar{v})$ denote the gradient of $f$ with respect to $u$ and $v$, respectively, computed at $(\bar{u},\bar{v})$.

2 Identification of Sparse AR Graphical Models

Assume to collect the data $y^{N}:=\{y(1),y(2)\ldots y(N)\}$ generated by the AR Gaussian discrete-time zero mean full rank stationary stochastic process $y=\{\,y(t),\;t\in\mathbb{Z}\}$ defined as

[TABLE]

$y(t)$ takes values in $\mathbb{R}^{m}$, $A_{k}\in\mathbb{R}^{m\times m}$ and $e(t)$ is white Gaussian noise with covariance matrix $R\succ 0$. Both $A_{k}$ and $R$ are unknown. The order of the AR process, i.e. $n$, is assumed to be known. We want to estimate $A_{k}$ and $R$ using $y^{N}$. It is well known that an equivalent description of $y$ is given by its PSD

[TABLE]

where $R_{k}=\mathbb{E}[y(t+k)y(t)^{T}]$, with $k\in\mathbb{Z}$, is the covariance lags sequence. Accordingly, the aforementioned problem is equivalent to estimate $\Phi$ from $y^{N}$. There are situations (e.g. $m$ is large) in which one is interested to estimate a PSD reflecting only the most important conditional dependence relations among the variables of $y$. Let $I\subset\mathcal{V}$ be an arbitrary index set. We denote as $\chi_{I}=\overline{\mathrm{span}}\{\,y_{j}(t)\hbox{ s.t. }\,j\in I\,,\,t\in\mathbb{Z}\,\}$ the closure of the vector space of all finite linear combinations (with real coefficients) of $y_{j}(t)$ with $j\in I$ and $t\in\mathbb{Z}$. Let $j\neq h$, we say that $y_{j}$ and $y_{h}$ are conditionally independent given the other variables if $\chi_{\{j\}}\,\bot\,\chi_{\{h\}}\,|\,\chi_{\mathcal{V}\setminus\{j,h\}}.$ These conditional dependence relations define an interaction graph $\mathcal{G}(\mathcal{V},\mathcal{E})$ where $\mathcal{V}$ and $\mathcal{E}$ denote the set of nodes and edges, respectively. More precisely, the nodes represent the variables $y_{1},y_{2}\ldots y_{m}$ and the lack of an edge means conditional independence (Brillinger, 1996):

[TABLE]

An example of graphical model is provided in Figure 1 (left).

Dahlhaus (2000) proved that $y_{j}$ and $y_{h}$ are conditionally independent if and only if $(\Phi^{-1}(e^{i\vartheta}))_{jh}=0$, $\forall\,\vartheta\in[-\pi,\pi].$ This characterization allows to infer conditional independence relations by promoting sparsity in the estimation of the inverse PSD for model (2). Since $y$ is an AR process of order $n$, we can parametrize its PSD as $\Phi=\Sigma^{-1}$ where $\Sigma=S_{0}+\frac{1}{2}\sum_{k=1}^{n}S_{k}e^{-ik\vartheta}+S_{k}^{T}e^{ik\vartheta}\in\mathcal{Q}_{m,n}$ and we define $S=[\,S_{0}\,S_{1}\,\ldots\,S_{n}\,]\in\mathbf{M}_{m,n}$. Then, a regularized maximum likelihood (ML) estimator of $\Sigma$, and thus of $\Phi$, is given by solving (Songsiri & Vandenberghe, 2010):

[TABLE]

The term $\ell(y^{N};\Sigma)$ is an approximation of the negative log-likelihood of $y(n+1)\ldots y(N)$ given $y(1)\ldots y(n)$ under model (2):

[TABLE]

where

[TABLE]

and $c$ is a term not depending on $y^{N}$ and $\Sigma$. It is worth noting that $\hat{R}_{k}$ represents an estimate of $R_{k}$ from data $y^{N}$, indeed $\Phi_{y^{N}}$ is the truncated periodogram of $\Phi$ computed from $y^{N}$. Under the assumption that $y$ is a full rank process, then for $N$ sufficiently large we have that $\mathrm{T}(\hat{R})\succ 0$ with high probability. Accordingly, throughout the paper we make the assumtion that $\mathrm{T}(\hat{R})\succ 0$. The penalty term $\bar{h}_{\infty}\,:\,\mathcal{Q}_{m,n}\longrightarrow\mathbb{R}$

[TABLE]

with

[TABLE]

encourages a common sparsity pattern (i.e. group sparsity) on the coefficients $S_{k}$ of $\Sigma$. Finally, $\gamma>0$ is the regularization parameter. Since $\Sigma\in\mathcal{Q}_{m,n}$ and in view of (1), Problem (2) can be rewritten in terms of a matrix $X\in\mathbf{Q}_{m(n+1)}$ (Songsiri & Vandenberghe, 2010):

[TABLE]

where $\ell(y^{N};X)=-\frac{N-n}{2}\log|X_{00}|+\frac{N-n}{2}\mathop{\rm tr}(\mathrm{T}(\hat{R})X)+c$; function $h_{\infty}\,:\,\mathbf{M}_{m,n}\longrightarrow\mathbb{R}$ is defined as (8) but now the domain is replaced by $\mathbf{M}_{m,n}$. Finally, the optimal solution $\hat{\Phi}=\hat{\Sigma}^{-1}$ is such that $\hat{\Sigma}=\Delta\hat{X}\Delta^{*}$. It is clear that $\hat{\Phi}^{-1}$ depends on $\gamma$. Songsiri & Vandenberghe (2010) proposed to select $\gamma$ by computing $j$ values of $\gamma$ according to the so called “trade-off” curve. Then the corresponding candidate models are ranked by using a BIC criterium. It is worth noting that the latter applies a thresholding on the partial coherence of the estimated PSD in order to measure the complexity (in terms of number of edges) of the candidate model.

3 Iterative Reweighted Method

In this section we investigate the possibility to modify the penalty term $\bar{h}_{\infty}$ in (2) in such a way to improve the ability to estimate the support of $\Phi^{-1}$. Notice that $h_{\infty}$ can be understood as the $\ell_{1}$ norm of the vector $[\,q_{11}(\Sigma)\,q_{12}(\Sigma)\ldots q_{nn}(\Sigma)\,]^{T}$ which represents the convex surrogate of the corresponding $\ell_{0}$ norm. As highlighted in Candès et al. (2008): “A key difference between the $\ell_{0}$ and the $\ell_{1}$ norms is the dependence on magnitude: larger coefficients are penalized more heavily in the $\ell_{1}$ norm than smaller coefficients, unlike the more democratic penalization of the $\ell_{0}$ norm”. We address this imbalance by considering a weighted penalty function $\bar{h}_{W}(\Sigma)=\sum_{j\geq h}\gamma_{jh}q_{jh}(\Sigma)$ where $\Sigma\in\mathcal{Q}_{m,n}$ and $\gamma_{jh}>0$. Note that, we penalize also the entries in the main diagonal of $\Sigma$. The latter regularization is not imposed in order to obtain sparsity in the main diagonal of $\Sigma$, otherwise constraint $\Sigma\succ 0$ is no longer satisfied, but rather it is imposed in order to reduce the variance of the estimator for the variables in the main diagonal of $\Sigma$. Accordingly, we will expect to find $\gamma_{jj}>0$ but smaller than $\gamma_{jh}$ with $j\neq h$. Thus, we modify Problem (2) as

[TABLE]

Following the same reasoning in Songsiri & Vandenberghe (2010) we can reformulate Problem (3) as follows:

[TABLE]

where $h_{W}\,:\,\mathbf{M}_{m,n}\rightarrow\mathbb{R}$ is defined as $\bar{h}_{W}$ but now the domain is $\mathbf{M}_{m,n}$. It is worth noting that (3) and (3) are equivalent provided that the optimal solution $\hat{X}$ is such that $\Delta\hat{X}\Delta^{*}\succ 0$.

Theorem 1**.**

The dual problem of (3) is

[TABLE]

The latter does admit a unique solution $\hat{X}$ such that $\Delta\hat{X}\Delta^{*}\succ 0$. Accordingly (3) and (3) are equivalent.

Theorem 1 is important not only for establishing the existence of the solution of (3) and (3), but it provides also a tractable formulation for the computation of the solution, see Section 6.1.

It remains to select a suitable set of weights $\gamma_{jh}$. The latter should counteract the magnitude imbalance which characterizes the $\ell_{1}$-like norm; more precisely, $\gamma_{jh}$ should be inversely proportional to $q_{jh}(\Sigma)$. For instance, we could take $\gamma_{jh}=(q_{jh}(\Sigma)+\varepsilon_{S})^{-1}$ with $\varepsilon_{S}>0$ sufficiently small. However, $\Sigma$ is unknown. Accordingly, we propose an iterative reweighted algorithm (Wipf & Nagarajan, 2010) which constructs iteratively a set of weights by using the information from the current estimate of $\Sigma$ (i.e. $X$), see Algorithm 1

where we recall that $S_{k}=\mathrm{D}_{k}(X)$. Parameter $\varepsilon_{S}>0$ in Step 6 ensures that a zero-valued entry in positions $(j,h)$ and $(h,j)$ does not strictly prohibit a nonnull estimate at the next step. Although Algorithm 1 has been introduced by an heuristic reasoning, it can be interpreted as the Majorization-Minimization (MM) algorithm for the problem

[TABLE]

The latter provides a regularized ML estimator of $\Sigma$ through the relation $S_{k}=\mathrm{D}_{k}(X)$ and $S_{k}$ are the coefficients of $\Sigma\in\mathcal{Q}_{m,n}$. The log-sum penalty induces sparsity on $\Sigma$. It is well known that such a penalty outperforms the corresponding $\ell_{1}$-like norm for estimating the correct sparsity pattern, see Candès et al. (2008). On the other hand, the log-sum penalty is concave making Problem (3) nonconvex. In order to see the aforementioned connection, we can rewrite (3) as

[TABLE]

Indeed, if $\hat{X}$ is the optimal solution of (3), then the optimal solution of (3) is $\hat{X}$ and $\hat{u}_{jh}=q_{jh}(\mathrm{D}(\hat{X}))$ with $j\geq h$. Let $u\in\mathbb{R}^{m(m+1)/2}$ be the vector obtained by stacking $u_{jh}$ with $j\geq h$. We define

[TABLE]

The latter is majorized by:

[TABLE]

More precisely, the first term in $f$, which is a concave function in $u$, is majorized by the tangent at $\hat{u}^{(l)}$. Then, the MM algorithm is

[TABLE]

Removing the terms in $g$ not depending on $u$ and $X$, we obtain

[TABLE]

Since we have $u_{jh}=q_{jh}(X)$ for the optimal solution of (3), then the latter problem is equivalent to

[TABLE]

Defining

[TABLE]

and substituting it in (3) we obtain (3) where $\gamma_{jh}$ has been replaced by $\hat{\gamma}_{jh}^{(l+1)}$. In other words, (3) is equivalent to: compute $\hat{\gamma}_{jh}^{(l+1)}$ as in (21) and then solve (3) with $\gamma_{jh}=\hat{\gamma}_{jh}^{(l+1)}$. The latter procedure coincides with the iterative reweighted scheme in Algorithm 1.

Proposition 2**.**

Consider the sequence $\hat{X}^{(l)}$ generated by Algorithm 1. Then, $\hat{X}^{(l)}$ converges to the set of stationary points of Problem (3).

Simulation evidence showed that the proposed iterative reweighted algorithm does not perform well, that is the estimated PSD is not close to the actual one, and neither the sparsity pattern of the inverse. Even if the updating rule for $\gamma_{jh}$s seems reasonable, it is not the best choice that we can apply. In the next section we will see how to find a better update.

4 A Bayesian Perspective

Till now the problem of estimating $\Sigma$ sparse has been considered according to the Fisherian perspective that is $\Sigma$ is an unknown but fixed function in $\mathcal{Q}_{m,n}$, i.e. the parameters $S_{k}$, $k=0\ldots n$, characterizing $\Sigma$ are unknown but fixed quantities. In this section we propose a method which is based on a Bayesian perspective that is $\Sigma$ is a stochastic process taking values in $\mathcal{Q}_{m,n}^{+}:=\{\Sigma\hbox{ s.t. }\Sigma\succ 0\}$. This means that the parameters characterizing $\Sigma$ are random variables with a suitable PDF or simply prior. Let $p(\Sigma)$ be the prior of $\Sigma$. We recall that $(\Sigma)_{jh}$ is the entry of $\Sigma$ in position $(j,h)$:

[TABLE]

We assume that $(\Sigma)_{jh}$s are independent each other, accordingly

[TABLE]

and

[TABLE]

where $c_{jh}$ is the normalizing constant and $\gamma_{jh}>0$, $j\geq h$, are referred to as hyperparameters. The latter are the parameters characterizing the prior. Therefore, we have

[TABLE]

where $\boldsymbol{\gamma}\in\mathbb{R}^{m(m+1)/2}$ denotes the hyperparameters vector containing $\gamma_{jh}$ with $j\geq h$. The negative log-likelihood of $y^{N}$ and $\Sigma$ takes the form:

[TABLE]

It is clear that the negative log-conditional PDF $-\log p(y^{N}|\Sigma)$ does coincide with (6), thus

[TABLE]

where the last term does not depend on $\Sigma$. We conclude that the MAP estimator of $\Sigma$ is given by (3). This result is not surprising, indeed it is well known that MAP estimators in a Bayesian perspective can be interpreted as regularized estimators in the Fisherian perspective. The substantial difference between the two perspectives is that the former provides the way to estimate the hyperparameters vector $\boldsymbol{\gamma}$. We define as the negative log-marginal likelihood:

[TABLE]

An estimate of $\boldsymbol{\gamma}$ is given by the empirical Bayes approach (Friedman et al., 2001):

[TABLE]

Then, the MAP estimator of $\Sigma$ is given by (3) with $\boldsymbol{\gamma}=\hat{\boldsymbol{\gamma}}$. However, it is not possible to find an analytical expression for $\ell(y^{N};\boldsymbol{\gamma})$ making challenging the optimization of $\boldsymbol{\gamma}$.

Proposition 3**.**

Consider the prior of $\Sigma$ defined in (22). Then, we have

[TABLE]

where the terms $\upsilon_{jh}$ in the relation above do not depend on $\boldsymbol{\gamma}$.

Accordingly,

[TABLE]

where $c$ is a term not depending on $y^{N}$, $\Sigma$ and $\boldsymbol{\gamma}$. An alternative simplified approach to estimate $\boldsymbol{\gamma}$ is the generalized maximum likelihood (GML) method, Zhou et al. (1997): instead of maximizing $\ell(y^{N};\boldsymbol{\gamma})$ with respect to $\boldsymbol{\gamma}$, the latter is computed with $\Sigma$ as the pair $(\hat{\boldsymbol{\gamma}},\hat{\Sigma})$ that jointly minimizes $\tilde{\ell}(y^{N},\Sigma;\boldsymbol{\gamma})$. Then, the optimization can be performed in a two-step algorithm:

[TABLE]

Step (31) is the MAP estimator of $\Sigma$ given the current choice of $\boldsymbol{\gamma}$ which is equivalent to (3). Step (32) is the estimator of $\boldsymbol{\gamma}$ using the current MAP estimate of $\Sigma$ as a direct observation. Notice that (32) is equivalent to

[TABLE]

where we recall that $\hat{S}^{(l)}=\mathrm{D}(\hat{X}^{(l)})$.

Proposition 4**.**

Under the assumption that $q_{jh}(\hat{S}^{(l)})>0$, it holds that

[TABLE]

To deal also with the case that $q_{jh}(\hat{S}^{(l)})=0$, we consider the modified updating

[TABLE]

with $\varepsilon_{S}>0$. It is not difficult to see that the latter modification is equivalent to assume that $\gamma_{jh}$s are modeled as independent random variables with exponential hyperprior: $p(\gamma_{jh})=\varepsilon_{S}e^{-\varepsilon_{S}\gamma_{jh}}$. We conclude that the GML approach is equivalent to Algorithm 1 wherein Step 6 is now replaced by (43). Finally, it is not difficult to see that the sequence $\hat{X}^{(l)}$ generated by the GML method converges to the set of stationary points of the Problem (3) where the log-sum penalty now is

[TABLE]

5 Identification of Latent-variable Graphical Models

Consider a Gaussian discrete-time zero mean full rank stationary stochastic process $z=\{\,z(t),\;t\in\mathbb{Z}\,\}$. $z$ is composed by $m$ manifest variables and $r$ latent variables, so that $z=[\,(y)^{T}\,(x)^{T}\,]^{T}$ with $y=[\,y_{1}\ldots y_{m}\,]^{T}$ and $x=[\,x_{1}\ldots x_{r}\,]^{T}$. We assume that $y$ is an AR process of order $n$ (known) and dimension $m$ (known). $r$ as well as the PSD of $z$ are unknown. We assume to collect the data $y^{N}:=\{\,y(1),y(2)\ldots$ $y(N)\,\}$ from the manifest process $y$. We want to estimate the PSD of $y$, say $\Phi$, in such a way that it corresponds to a latent-variable graphical model. The latter is an interaction graph with two layers: latent nodes $x_{1}\ldots x_{r}$ are in the upper level, manifest nodes $y_{1}\ldots y_{m}$ are in the lower level, $r\ll m$ and there are few edges among the manifest nodes. An example of a latent-variable graphical model is provided in Figure 1 (right). The powerfulness of latent-variable graphical models is that the introduction of few latent variables (in respect to the manifest ones) may reduces drastically the conditional dependences among the manifest variables. Accordingly, in such graphs we expect that the interdependence relations among the manifest variables are mainly explained by few and common latent variables. In Zorzi & Sepulchre (2016) it has been shown that $\Phi$ corresponds to a latent-variable graphical model if it admits the following decomposition

[TABLE]

where

[TABLE]

$\Sigma\succ 0$ is sparse and its support reflects the conditional dependence relations among the manifest variables. $\Lambda\succeq 0$ is low rank and its rank is equal to $r$, i.e the number of latent variables. The approximate negative log-likelihood of $y(n+1)\ldots y(N)$ given $y(1)\ldots y(n)$ under the aforementioned model is:

[TABLE]

where $\ell(y^{N};\cdot)$ has been defined as in (6). Accordingly, the regularized ML estimator of $\Sigma$ and $\Lambda$ is (Zorzi & Sepulchre, 2016):

[TABLE]

Here, $\bar{r}_{\star}(\Lambda)=\mathop{\rm tr}\int\Lambda$ is the nuclear norm of $\Lambda\in\mathcal{Q}_{m,n}$ with $\Lambda\succeq 0$. $\gamma_{S},\gamma_{L}>0$ are the regularization parameters for $\Sigma$ and $\Lambda$, respectively. The term $\bar{h}_{\infty}(\Sigma)$, defined in (8), induces sparsity on $\Sigma$. The term $\bar{r}_{\star}(\Lambda)$ induces low rank on $\Lambda$, indeed it has been shown that such a function is the convex envelop of $\mathrm{rank}(\Lambda)$, see Proposition 3.1 in Zorzi & Sepulchre (2016). It is worth noting that the estimate of $r$ is given by the numerical rank of $\hat{\Lambda}$. Since $\Sigma-\Lambda$ and $\Lambda$ belong to $\mathcal{Q}_{m,n}$ and in view of (1), we can define $X,H\in\mathbf{Q}_{m(n+1)}$ such that

[TABLE]

and $X,H\succeq 0$. It is possible to prove that Problem (5) can be rewritten in terms of $X$ and $H$ as follows:

[TABLE]

where $\ell(y^{N};X)$ does coincide with the first term in the objective function of (2). Also in this case the optimal solution $(\hat{\Sigma},\hat{\Lambda})$ depends on the regularization parameters $\gamma_{S}$, $\gamma_{L}$. The values $(\gamma_{S},\gamma_{L})$ can be selected by considering a 2-dimensional grid. Then, the candidate models can be ranked by using a BIC criterium similar to the one in Songsiri & Vandenberghe (2010). The latter measures the complexity by thresholding the partial coherence of the sparse component and the singular values of the low rank component. Alternatively, Zorzi & Sepulchre (2016) proposes a score function based on the Kullback-Leibler divergence to rank the candidate models.

Similarly to the case without latent variables, the weaknesses of the regularization Problem (5) are that: (i) the nonnull entries in $\Sigma$ are penalized in a different way; (ii) the nonnull eigenvalues of $\Lambda$, understood as functions over the unit circle, are penalized in a different way. The sparse regularizer can be made “democratic” in the same way of before. Regarding the low rank part, we can replace $\bar{r}_{\star}(\Lambda)$ with the “weighted” penalty:

[TABLE]

where the weight matrix $Q\in\mathbf{Q}_{m}$ is such that $Q\succ 0$. Thus, we consider the problem:

[TABLE]

Using the parametrization in (48), we have that

[TABLE]

where we exploited the well known identity $\int e^{ik\theta}=\delta_{k}$ and $\delta_{k}$ denotes the Kronecker delta function. Thus, we consider the problem:

[TABLE]

Note that, (5) and (5) are equivalent provided that $\Delta\hat{X}\Delta^{*}\succ 0$.

Theorem 5**.**

The dual problem of (5) is

[TABLE]

The latter does admit solution $(\hat{X},\hat{H})$ such that $\hat{X}$ is unique and $\Delta\hat{X}\Delta^{*}\succ 0$. Accordingly (5) and (5) are equivalent.

Theorem 5 is important not only for establishing the existence of the solution of (5) and (5), but it provides also a tractable formulation for the computation of the solution, see Section 6.1. It is worth noting that Problem (5) may have more than one solution. On the other hand, if we compute an optimal solution of (5), then from the latter we can recover the solution of (5) by solving a system of linear equations in $H$ (a similar idea has been used in Section III.C in Zorzi & Sepulchre (2016)). The uniqueness of the solution of this system of linear equations is guaranteed provided that: (i) $Q$ has a sufficient number of eigenvalues which are sufficiently large; (ii) there is a sufficient number of $\gamma_{jh}$s with $j\neq h$ taking sufficiently large values. The latter implies the uniqueness of the solution to (5) and thus the uniqueness of the one to (5).

To select a suitable set of $\gamma_{jh}$ and a suitable $Q$, we design a reweighted algorithm which exploits a Bayesian perspective. We model $\Sigma$ and $\Lambda$ as stochastic processes taking values in $\mathcal{Q}_{m,n}^{+}$ and in the closure of $\mathcal{Q}_{m,n}^{+}$, respectively. Let $p(\Sigma)$ and $p(\Lambda)$ denote the prior of $\Sigma$ and $\Lambda$, respectively. We assume that $\Sigma$ and $\Lambda$ are independent, that is the joint PDF of $\Sigma$ and $\Lambda$ is such that $p(\Sigma,\Lambda)=p(\Sigma)p(\Lambda)$. We set the prior of $\Sigma$ as in (22). Regarding $\Lambda$, $L_{k}$s are modeled as independent random matrices. More precisely, we model $L_{0}$ as a Wishart random matrix with $m+1$ degrees of freedom and variance $mQ^{-1}$:

[TABLE]

where $Q\in\mathbf{Q}_{m}$ such that $Q\succ 0$ and $G_{m}$ is the multivariate gamma function. Finally, we attach an uninformative prior on $L_{k}$ with $k\geq 1$. Then, the negative log-likelihood of $y^{N}$, $\Sigma$ and $\Lambda$ is defined as

[TABLE]

where $c$ is a constant term not depending on $y^{N}$, $\Sigma$, $\Lambda$, $\boldsymbol{\gamma}$ and $Q$. Note that, $-\log p(y^{N}|\Sigma,\Lambda)$ does coincide with $\ell(y^{N};\Sigma,\Lambda)$ defined in (46). Also in this case, it is not possible to find an analytical expression for the negative log-marginal likelihood. Therefore, we have the following upper bound for $\ell(y^{N},\Sigma,\Lambda;\boldsymbol{\gamma},Q)$:

[TABLE]

Then, according to the GML approach we consider the following two-step algorithm:

[TABLE]

Clearly, (59) is the MAP estimator of $\Sigma$ and $\Lambda$ given the current choice of $\boldsymbol{\gamma}$ and $Q$, while (60) is the estimator of $\boldsymbol{\gamma}$ and $Q$ using the current MAP estimate of $\Sigma$ and $\Lambda$ as a direct observation. The objective function in (60) can be split in two terms depending on $\boldsymbol{\gamma}$ and $Q$, respectively. Accordingly, $\hat{\boldsymbol{\gamma}}^{(l+1)}$ takes a form similar to the one in (40). Regarding $Q$, we have:

[TABLE]

Proposition 6**.**

Under the assumption that $\hat{L}_{0}^{(l)}\succ 0$, we have that

[TABLE]

To deal also with the case $\hat{L}_{0}^{(l)}$ singular matrix we consider the modified updating rule:

[TABLE]

where $\varepsilon_{L}>0$. It is not difficult to see that the latter modification is equivalent to assume that $Q$ is modeled as a Wishart random matrix with $m+1$ degrees of freedom and variance $m\varepsilon_{L}^{-1}I$, i.e. $p(Q)=e^{-\varepsilon_{L}\mathop{\rm tr}(Q)}/(\sqrt{2^{(m+1)m}\varepsilon_{L}^{-m(m+1)}}G_{m}((m+1)/2))$, which is assumed to be independent of $\boldsymbol{\gamma}$. Algorithm 2 describes the corresponding procedure which is clearly an iterative reweighted algorithm. Here, we recall that $\hat{S}^{(l)}=\mathrm{D}(\hat{X}^{(l)}+\hat{H}^{(l)})$ and $\hat{L}_{0}^{(l)}=\mathrm{D}_{0}(\hat{H}^{(l)})$.

It is not difficult to see that Algorithm 2 can be interpreted as the MM algorithm for solving the regularized ML problem:

[TABLE]

where

[TABLE]

Beside the sparse-inducing part (already analyzed in Section 4), it is known that the log-det penalty outperforms the nuclear norm for estimating the correct rank, see Fazel et al. (2003).

Corollary 7**.**

Consider the sequence $(\hat{X}^{(l)},\hat{H}^{(l)})$ generated by Algorithm 2 and assume that the solution of (5), with weights $\hat{\boldsymbol{\gamma}}^{(l)}$ and $\hat{Q}^{(l)}$, is always unique. Then, $(\hat{X}^{(l)},\hat{H}^{(l)})$ converges to the set of stationary points of Problem (5).

6 Simulation Results

In this section we test the performance of the proposed empirical Bayes estimators for learning AR graphical models.

6.1 Implementation details

The empirical Bayes estimate of a sparse AR graphical model is obtained by running Algorithm 1 wherein Step 6 is replaced by the updating rule (43). The most challenging part is the computation of $\hat{X}^{(l)}$ (Step 4) which is obtained through the dual problem (1). The optimal solution of the latter is computed by implementing a first-order projection algorithm. The empirical Bayes estimate of a latent-variable AR graphical model is obtained by running Algorithm 2. Here, the most challenging part is the computation of $(\hat{X}^{(l)},\hat{H}^{(l)})$, in Step 4, which is obtained through the dual problem (5). The optimal solution of the latter is computed by using the alternating direction method of multipliers (ADMM). More precisely, we introduce an auxiliary variable $K=I\otimes Q+\mathrm{T}(Z)$. In the first step, we minimize the augmented Lagrangian with respect to $W$ and $Z$ under the first four constraints of (5). This minimization is performed by using a first-order projection algorithm. In the second step, we minimize the augmented Lagrangian with respect to $K$ under the last constraint of (5). This optimal solution can be computed in closed form. Finally, the corresponding Matlab functions are available at https://github.com/MattiaZ85/EBL-AR-GM.

6.2 Identification of sparse graphical models

We consider two Monte Carlo experiments which are structured as follows. We generate 200 AR models of dimension $m=30$, of order $n$ and whose fraction of nonnull entries in the inverse PSD is equal to $0.1$. The position of such nonnull entries is chosen randomly for each model. For each model we generate a finite data sequence of length $N$ and we consider the following estimators:

• •

TD9 is the estimator proposed in Songsiri & Vandenberghe (2010) where the number of points for tracing the trade-off curve is set equal to $9$ and the threshold for the partial coherence is set equal to $0.1$;

• •

TD17 is the estimator proposed in Songsiri & Vandenberghe (2010) where the number of points for tracing the trade-off curve is set equal to $17$ and the threshold for the partial coherence is set equal to $0.1$;

• •

RW is the empirical Bayes estimator that we have proposed in Section 4 where $\varepsilon_{S}=10^{-3}$. The hyperparameters vector $\boldsymbol{\gamma}$ is initialized with (43) where $\hat{S}^{(l)}=\mathrm{D}(X_{B})$, $X_{B}\in\mathbf{Q}_{m(n+1)}$ is such that $\hat{\Phi}_{B}=(\Delta X_{B}\Delta^{*})^{-1}$ and $\hat{\Phi}_{B}$ is the Burg estimator computed from the data111It is worth recalling that the Burg method provides an estimate of the PSD which corresponds to an AR model of a certain order, fixed by the user. In our case, the latter is fixed equal to $n$; in this way $\hat{\Phi}_{B}=(\Delta X_{B}\Delta^{*})^{-1}$.. The threshold for the stopping condition is set equal to $\varepsilon=10^{-4}$ and the maximum number of iterations is set equal to $l_{\mathrm{MAX}}=50$.

For each estimator:

• •

we compute the relative error in the estimation of the inverse PSD coefficients as $e=\|\hat{S}-S\|^{2}/\|S\|^{2}$ where $\hat{S}=[\,\hat{S}_{0}\,\hat{S}_{1}\,\ldots\,\hat{S}_{n}]\in\mathbf{M}_{m,n}$ are the coefficients of the inverse of $\hat{\Phi}$ (estimated PSD), while $S=[\,S_{0}\,S_{1}\,\ldots\,S_{n}]\in\mathbf{M}_{m,n}$ are the ones of the inverse of $\Phi$ (true PSD);

• •

we compute the fraction of null and nonnull misplaced entries in the inverse of the estimated PSD with respect to the inverse of the true PSD as $e_{SP}=\#\hat{E}_{m}/(m(m-1)/2)$ where $\#\hat{E}_{m}$ denotes the total number of null and nonnull misplaced entries in the inverse of the estimated PSD.

In the first experiment, the AR models are of order $n=1$ and we have considered two different lengths of the data, that is $N=500$ and $N=1000$. The boxplots of $e$ and $e_{SP}$ are depicted in Figure 2. As we can see, RW outperforms TD9 and TD17 both in terms of estimation error and of misplaces entries. It is worth noting that TD9 and TD17 perform in the same way: this confirms the observation that for tracing accurately the trade-off curve it is required just a small number of points (Songsiri & Vandenberghe, 2010). We have noticed that TD9 and TD17 estimate the inverse PSD too sparse in respect to the true one. For this reason, we also have considered TD9 and TD17 with the threshold for the partial coherence equal to $0.05$: the resulting performances are worse than the ones with $0.1$. In the second experiment, the AR models are of order $n=2$ and we have considered the two different lengths of the data as before. The boxplots of $e$ and $e_{SP}$ are depicted in Figure 3. The results are similar to the ones of the previous experiment: RW outperforms the other two estimators; TD9 and TD17 perform in the same way.

6.3 Identification of latent-variable AR graphical models

We consider two Monte Carlo experiments which are structured as follows. We generate 200 AR models of dimension $m=30$ of order $n=2$ whose inverse PSD has a sparse plus low rank decomposition, i.e. $\Phi^{-1}=\Sigma-\Lambda$ with $\Sigma$ sparse and $\Lambda$ low rank. The fraction of nonnull entries in the sparse part is denoted by $s$; the position of the nonnull entries is chosen randomly for each model. The rank of $\Lambda$ is denoted by $r$; for each model, the image of $\Lambda$ is chosen randomly in such a way that $\Sigma-\Lambda$ is positive definite over the unit circle. For each model we generate a finite data sequence of length $N=1000$ and we consider the following estimators:

• •

TD-KL is the estimator proposed in Zorzi & Sepulchre (2016) where the model is chosen using the score function based on the Kullback-Leibler divergence; the threshold for the partial coherence (sparse component) is set equal to $0.1$; the threshold for the normalized singular values (low rank component) is set equal to $0.1$; the regularization grid (16 points) is chosen in such a way that the candidate model, say $\hat{\Phi}=(\hat{S}-\hat{L})^{-1}$, ranges from $\hat{S}$ diagonal - $\mathrm{rank}(\hat{L})$ large to $\hat{S}$ moderately sparse and $\mathrm{rank}(\hat{L})=0$.

• •

TD-BIC is the estimator proposed in Zorzi & Sepulchre (2016) where the model is chosen using the BIC score function; the threshold parameters and the regularization grid are chosen as in the previous estimator.

• •

RW is the empirical Bayes estimator that we have proposed in Section 5 where $\varepsilon_{S}=10^{-3}$ and $\varepsilon_{L}=10^{-3}$. The hyperparameters are initialized as follows. Let $\hat{\Phi}_{B}=(\Delta X_{B}\Delta^{*})^{-1}$ be the Burg estimator of the PSD computed from the data. Then, the initial value for $\boldsymbol{\gamma}$ is given by (43) where $\hat{S}^{(l)}=(1+\alpha)\mathrm{D}(X_{B})$, the initial value for $Q$ is given by (63) where $\hat{L}^{(l)}=\alpha\mathrm{D}_{0}(X_{B})$ and $\alpha=0.1$. The threshold for the stopping condition is set equal to $\varepsilon=10^{-4}$ and the maximum number of iterations is set equal to $l_{\mathrm{MAX}}=50$.

For each estimator:

• •

we compute the relative error in the estimation of the inverse PSD coefficients as $e=\|\hat{S}-\hat{L}-(S-L)\|^{2}/\|S-L\|^{2}$ where $\hat{S}=[\,\hat{S}_{0}\,\hat{S}_{1}\,\ldots\,\hat{S}_{n}]\in\mathbf{M}_{m,n}$ are the coefficients of $\hat{\Sigma}$, while $\hat{L}=[\,\hat{L}_{0}\,\hat{L}_{1}\,\ldots\,\hat{L}_{n}]\in\mathbf{M}_{m,n}$ are the ones of $\hat{\Lambda}$. $\hat{\Phi}=(\hat{\Sigma}-\hat{\Lambda})^{-1}$ denotes the estimated PSD. In a similar way, $S$ and $L$ contain the coefficients of $\Sigma$ and $\Lambda$ such that $\Phi=(\Sigma-\Lambda)^{-1}$ is the true PSD.

• •

we compute the relative error in the estimation of the sparse and low rank coefficients in the inverse PSD as $e_{SL}=(\|\hat{S}-S\|^{2}+\|\hat{L}-L\|^{2})/\|S-L\|^{2}$;

• •

we compute the complexity of the estimated model as $C=\#p/(m(m+1)/2+m^{2}n)$ where $\#p$ is the number of parameters needed to characterize the model $\hat{\Phi}=(\hat{\Sigma}-\hat{\Lambda})^{-1}$. For instance if $\hat{\Sigma}$ has $s$ nonnull entries and $\mathrm{rank}(\Lambda)=r$, then $\#p=(s+m)/2+sn+rm(n+1)$. The denominator in $C$ is the number of coefficients needed in an unstructured model, i.e. the corresponding interaction graph is complete.

In the first experiment, the AR models are with $s=0.1$ and $r=2$. The relative errors $e$ and $e_{SL}$ are depicted in Figure 4. As we can see, RW provides the best performance while TD-KL the worst performance. The average complexity of the estimated models for each estimator (and denoted by $\bar{C}$) is shown in the table on the left of Figure 6. As we can see, the average complexity of the models estimated by RW is closer to the true one than the one of TD-KL and TD-BIC. More precisely, TD-KL and TD-BIC estimate models which are extremely simple, in respect to the true complexity, and the price to pay is the inferiority in respect to $e$ and $e_{SL}$. In the second experiment, the AR models are with $s=0.1$ and $r=5$. The relative errors $e$ and $e_{SL}$ are depicted in Figure 5, while the average complexity of the estimated models for each estimator is shown in the table on the right of Figure 6. As we can see, the conclusions of the previous case still hold. Finally, we have tested RW for different values of $\alpha$ such that $0.08\leq\alpha\leq 0.15$. We have noticed that RW provides similar performances of before.

7 Conclusion

We have analyzed the problem of estimating sparse and latent-variable AR graphical models. These two problems are traditionally solved by using $\ell_{1}$-like and nuclear norm-like regularizes. The latter introduce a magnitude imbalance in the optimization problem which produces an estimator with a remarkable mean squared error. We have proposed two empirical Bayes estimators which counteract such an imbalance. The hyperparameters of these estimators are computed by the generalized maximum likelihood approach which leads to an iterative reweighted method. Simulation results showed the benefit to introduce a prior for the model that we have to estimate.

Appendix

Proof of Theorem 1

First, we can rewrite (3) by adding a new variable $Y\in\mathbf{M}_{m,n}$:

[TABLE]

The Lagrangian is

[TABLE]

where $U\in\mathbf{Q}_{m(n+1)}$ is such that $U\succeq 0$ and $Z\in\mathbf{M}_{m,n}$. We proceed to perform the unconstrained minimization of $\mathcal{L}$ with respect to the primal variables. It is not difficult to see that:

[TABLE]

Accordingly, the latter admits minimum if and only if

[TABLE]

moreover, in such a situation it takes value equal to zero. Thus,

[TABLE]

if (Proof of Theorem 1) holds, otherwise it takes $-\infty$. Regarding the minimization of $\mathcal{L}$ with respect to $X$: the minimum with respect to $X_{00}$ is $X_{00}=(\mathrm{T}(\hat{R})+\mathrm{T}(Z)-U)_{00}^{-1}$; regarding the other blocks $X_{jh}$, it is not difficult to see that $\mathcal{L}$ is bounded below if and only if

[TABLE]

We conclude that:

[TABLE]

if (Proof of Theorem 1)-(70) hold, otherwise it takes $-\infty$. Substituting $U$ with the new variable $W:=(\mathrm{T}(\hat{R})+\mathrm{T}(Z)-U)_{00}$, we obtain (1). Regarding the existence of the solution, notice that Problem (1) consists in maximizing $\log|W|$ subject to the constraints $W\succ 0$, \mathrm{T}(\hat{R})+\mathrm{T}(Z)\succeq\left[\begin{array}[]{cc}W&0\\ 0&0\end{array}\right], and $Z\in\mathbf{C}$ where $\mathbf{C}:=\{\,Z\in\mathbf{M}_{m,n}\hbox{ s.t. }(\ref{cond1})\hbox{ holds}\,\}$ which is a closed convex subset of $\{Z\in\mathbf{M}_{m,n}\hbox{ s.t. }\mathop{\rm tr}(Z_{0})$ $\hbox{ bounded}\}$. Then, the existence of the solution is guaranteed by applying a modified version of Lemma A.1 in Zorzi & Sepulchre (2016). More precisely, it is not difficult to see that the conclusion of such Lemma does not change by replacing the condition $\mathop{\rm tr}(Z_{0})=0$ with $\mathop{\rm tr}(Z_{0})$ bounded. The remaining part of the proof follows the same lines of the one of Proposition 3.3 in Zorzi & Sepulchre (2016). ∎

Proof of Proposition 2 and Corollary 7

We need the following result.

Lemma 7.1**.**

(Razaviyayn et al., 2013, Corollary 1)* Consider the problem*

[TABLE]

which is assumed to be feasible, i.e. there exists $\bar{u}\in\mathcal{C}$ such that $f(\bar{u})$ takes finite value, and $\mathcal{C}$ is a convex set. Let $\mathcal{M}$ denote the set of stationary points for (71). Consider the corresponding MM algorithm

[TABLE]

where $g(u;v)$ is such that $g(u;u^{(l)})\geq f(u^{(l)})$, $g(u^{(l)};u^{(l)})=f(u^{(l)})$, $\nabla_{u}g(u;u)=\nabla_{u}f(u)$ and $g(u;v)$ is continuous in $(u,v)$. If the level set $\mathcal{S}:=\{u\in\mathcal{C}\hbox{ s.t. }f(u)\leq f(\bar{u})\}$ is compact, then the sequence $u^{(l)}$ converges to the set $\mathcal{M}$, i.e. $\lim_{l\rightarrow\infty}\,d(x^{(l)},\mathcal{M})=0$ where $d(x^{(l)},\mathcal{M}):=\inf_{s\in\mathcal{M}}\|x-s\|$.

We proceed to prove Proposition 2 by applying Lemma 7.1 to Problem (3) and its MM algorithm (3). It is not difficult to see all the hypotheses on the approximating function $g$ are satisfied. It remains to prove the feasibility and compactness of the level set. Notice that (3) and (3) are equivalent, thus we can consider Problem (3). In our case, $\mathcal{C}=\{\,X\,\hbox{ s.t. }\;X_{00}\succ 0,\;X\succeq 0\,\}$ which is a convex set, and $f(X)=\sum_{j\geq h}\log(q_{jh}(\mathrm{D}(X))+\varepsilon_{S})+\ell(y^{N};X)$. Problem (3) is feasible indeed it is sufficient to pick $\bar{X}=I$. Consider the level set $\mathcal{S}:=\{\,X\in\mathcal{C}\;\hbox{ s.t. }\;\;f(X)\leq f(\bar{X})\,\}$. We show that $\mathcal{S}$ is closed and bounded, since it is a subset of a finite dimensional vector space, then $\mathcal{S}$ is compact. Consider a convergent sequence $X^{(l)}\in\mathcal{C}$, $l\in\mathbb{N}$, such that $X_{00}^{(l)}$ converges to a singular matrix, then $-\log|X_{00}^{(l)}|$ tends to infinity. Then, $f(X^{(l)})\geq m(m+1)/2\log\varepsilon_{S}-\log|X_{00}^{(l)}|\rightarrow\infty$. Accordingly, such a sequence cannot belong to $\mathcal{S}$, and thus $\mathcal{S}$ is a closed set. We proceed to show that unbounded sequences cannot belong to $\mathcal{S}$, i.e. $\mathcal{S}$ is bounded. Consider a convergent sequence $X^{(l)}\in\mathcal{C}$, $l\in\mathbb{N}$ such that $\|X^{(l)}\|\rightarrow\infty$. This means that $X^{(l)}$ has at least one eigenvalues tending to infinity as $l\rightarrow\infty$. In $f(X^{(l)})$ the term $\mathop{\rm tr}(\mathrm{T}(\hat{R})X^{(l)})\rightarrow\infty$ because $\mathrm{T}(\hat{R})\succ 0$. Moreover, it dominates the other two terms. We conclude that $f(X^{(l)})\rightarrow\infty$, and thus the sequence cannot belong to $\mathcal{S}$. We conclude that $\mathcal{S}$ is bounded, and thus compact. Thus, we can apply Corollary 1 and hence the claim is proved. Finally, Corollary 7 can be proved in a similar way. ∎

Proof of Proposition 3

Before to prove the statement we need the following Lemma.

Lemma 7.2**.**

Let $\gamma>0$. Consider the integral:

[TABLE]

Then, $\gamma^{-n}\leq\mathcal{I}\leq n^{n}\gamma^{-n}.$

Proof 7.3**.**

Notice that:

[TABLE]

accordingly

[TABLE]

Taking the integral with respect to $x_{1},\ldots,x_{k}$ over $[0,\infty)\times\ldots\times[0,\infty)$ in the above inequalities, we obtain a lower and an upper bound for $\mathcal{I}$, respectively:

[TABLE]

which concludes the proof. ∎We proceed to prove Proposition 3. Since $c_{jh}$ is the normalizing constant of $p_{\gamma_{jh}}((\Sigma)_{jh})$, then we have that

[TABLE]

where

[TABLE]

and

[TABLE]

In the case that $j\neq h$, $\tilde{c}_{jh}$ is an integral taken with respect to $2n+1$ variables. Since $\tilde{q}_{jh}(S)=\tilde{q}_{jh}(-S)$, then

[TABLE]

Hence, by Lemma 7.2 we have that

[TABLE]

In the case that $j=h$, we have that $(S_{k})_{jh}=(S_{k})_{hj}$ for $k\geq 0$, accordingly

[TABLE]

and the integration is taken with respect to $n+1$ variables. Along the same reasoning of before we have that

[TABLE]

In view of (73) and (74), we get the claim where $v_{jh}=2^{2n+1}(2n+1)^{2n+1}$ or $v_{jh}=2^{n+1}(n+1)^{n+1}$. ∎

Proof of Proposition 4

Consider the case $j=h$, we have to minimize the objective function $f_{jh}(\gamma_{jh}):=\gamma_{jh}q_{jh}(\hat{S}^{(l)})-(n+1)\log\gamma_{jh}$ under the constraint that $\gamma_{jh}>0$. Notice that $f_{jh}(\gamma_{jh})$ is strictly convex. Thus, the minimum point is given by setting equal to zero the first derivate: $q_{jh}(\hat{S}^{(l)})-(n+1)\gamma_{jh}^{-1}=0$ which is satisfied with $\gamma_{jh}=(n+1)/q_{jh}(\hat{S}^{(l)})$ and the constraint $\gamma_{jh}>0$ is satisfied. The proof for the case $j>h$ is similar. ∎

Proof of Theorem 5

The proof is similar to the one of Theorem 1. ∎

Proof of Proposition 6

The objective function in (61) is strictly convex in $Q$. Accordingly, the point of minimum is given by setting equal to zero the first variation of the objective function for any $\delta Q\in\mathbf{Q}_{m}$: $\mathop{\rm tr}(\hat{L}_{0}^{(l)}\delta Q-\frac{m+1}{2}Q^{-1}\delta Q)=0$, which implies that $\hat{L}_{0}^{(l)}-\frac{m+1}{2}Q^{-1}=0$. Therefore, under the assumption that $\hat{L}_{0}^{(l)}\succ 0$ we have (62).∎