Break recovery in graphical networks with D-trace loss

Much attention has been devoted to the development of methods for the detection of multiple change-points in the underlying data generating process of some realized sample. The model parameters are typically assumed to be constant in each regime but can experience “jumps” or “breaks” over time. The extraction of these breaks is usually performed through the application of some filtering techniques. To do so, a popular method is the fused LASSO which screens the existence of parameter breaks over the sample of observations. In particular, [17] developed a framework to recover multiple change-points for one-dimensional piece-wise constant signals. Then [4] extended this procedure to grouped parameters with the Group Fused LASSO in the context of autoregressive time series models. In the same spirit, [30] applied the Group Fused LASSO to linear regression models. This filtering technique has also been popular among the literature on change-points in panel data models: e.g., [31] aimed to identify structural changes in linear panel data regressions based on the adaptive Group Fused LASSO; using a similar filtering technique, [21] considered the estimation of structural breaks in panel models with interactive fixed effects.

This paper considers the detection of structural breaks in time series network, whose structure is represented by the corresponding precision matrix. Some works have been dedicated to change-point detection for precision matrix. [46] assumed that the covariance matrix evolves smoothly over time. Contrary to the framework of [20], which focused on the detection of multiple breaks at a node level through the Group Fused LASSO, [33] restricted their framework to a single change point which impacts the global network structure. Our aim is to detect multiple change-points that affect the whole network. Thus, our viewpoint is similar to [16] and [13]: these works considered the Group Fused Graphical LASSO (GFGL) to detect breaks in the precision matrix. The latter filtering technique is a mixture of the Group Fused LASSO regularization of [2] and the LASSO regularization applied to the Gaussian likelihood function. We propose an alternative estimator to the commonly employed penalized Gaussian likelihood, which builds upon the D-trace loss of [44]. The formulation of the latter is much simpler than the Gaussian likelihood (in particular, the absence of log-determinant), thus allowing for a more direct theoretical analysis and implementation procedure. The D-trace loss and its extensions have been applied to diverse problems: while [40] considered network change analysis between two samples of vector-valued observations, [19] adapted the latter framework to differential network analysis for fMRI data; in the context of compositional data, [41] and [43] considered modified versions of the D-trace loss to estimate the high-dimensional compositional precision matrix under sparsity constraint; using the matching of the coefficients of SVAR processes and the entries of the precision matrix of the observed sample, [29] considered the sparse estimation of the latter based on the D-trace loss. In addition to the detection of structural breaks in the network, we allow for sparse dependence in the entries of the precision matrix in each regime. This motivates the use of the LASSO penalization function applied to each coefficient of the piece-wise constant precision matrix.

The corresponding estimation problem formulated in (2.1), Section 2, with the D-trace loss, includes two tuning parameters. In practice, the optimal parameters are unknown a priori and are expected to be estimated upon solving the estimation problems with a series of tuning parameters. However, we show that the problem with general tuning parameters may be unbounded from below (and hence, does not have solutions). Consequently, to identify the optimal parameters, one may encounter estimation problems that are unbounded from below, and the unboundedness can be numerically difficult to detect. To address this issue, we introduce a new regularizer, thereby constructing a revised problem that consistently has solutions regardless of the choice of the tuning parameters. In addition to the existence of a solution, this revised problem also enjoys several desirable properties. On the one hand, if the solutions to the revised problem possess some easy-to-detect patterns, then the original problem may not have solutions, and we can further update the tuning parameters towards obtaining reliable estimators. On the other hand, if the patterns are not detected, then the solutions to the revised problem also solve the original problem. We adapt the celebrated alternating direction method of multipliers (ADMM) to solve the revised problem, with its convergence guaranteed by [10]. ADMM is a widely used algorithm for solving convex optimization problems with separable objectives and linear coupling constraints. Its classical version was first introduced by [15] and [12], with the convergence established in [11]. Then [7] showed that ADMM is equivalent to the proximal point algorithm applied to a certain maximal monotone operator. This insight led to the development of the first proximal ADMM by [6]. Building upon Eckstein’s work, [18] extended the proximal ADMM to include a broader class of proximal terms. This approach was further advanced by [10], which generalized the method to allow the use of semi-proximal terms, enhancing the algorithm’s applicability. More recently, considerable efforts have been made to further accelerate variants of ADMM: see, e.g., the Schur complement-based semi-proximal ADMM in [23]; the inexact symmetric Gauss-Seidel-based ADMM in [5], [24], [25], [39]; the accelerated linearized ADMM in [28], [37], [22]; the preconditioned ADMM in [36], [34], [35]; the accelerated proximal ADMM based on Halpern iteration in [42], [38].

Our contributions can be summarized as follows: we propose a new estimator for break detection in the precision matrix, the Group Fused D-trace LASSO (GFDtL); we derive the conditions for which all the break points and the precision matrices can be consistently estimated when the estimated breaks coincide with the true number of breaks; we provide a modified regularizer to ensure the existence of solutions to the revised problem, and show that these solutions either exhibit some easily verifiable patterns indicating the possible unsolvability of the original problem so that we can further update the tuning parameters towards obtaining reliable estimators, or solve the original problem if the patterns are not detected; an ADMM is adapted for implementation with convergence guarantees. The relevance of the novel estimator for change-points in time-varying networks compared to the standard Gaussian-based GFGL estimator is illustrated through simulations and real data experiments.

The rest of the paper is organized as follows. Section 2 details the framework and estimation procedure. Section 3 contains the asymptotic properties. Section 4 is devoted to the optimization aspect of the estimation procedure. Section 5 provides the implementation details. Section 6 contains some simulation experiments, a sensitivity and computational complexity analysis. A real data experiment is provided in Section 7. All the proofs of the main text and auxiliary results are relegated to the Appendices.

Notation: Throughout this paper, we use $V_{k}$ and $A_{kl}$ to denote the $k$-th element of a vector $V\in{\mathbb{R}}^{d}$ and the $(k,l)$-th element of a matrix $A\in{\mathbb{R}}^{m\times n}$, respectively. We write $A^{\top}$ (resp. $V^{\top}$) to denote the transpose of the matrix $A$ (resp. the vector $V$). For any square matrix $A\in{\mathbb{R}}^{n\times n}$, we write $\operatorname*{\mathrm{Sym}}(A)\coloneqq(A+A^{\top})/2$ to denote the symmetrization of $A$. The $p$-th order identity matrix is denoted by $I_{p}$. We denote by $\textbf{0}_{k\times l}\in{\mathbb{R}}^{k\times l}$ (resp. $\mathbf{1}_{k\times l}\in{\mathbb{R}}^{k\times l}$) the $k\times l$ zero matrix (resp. $k\times l$ matrix of ones). We write $\text{vec}(A)$ to denote the vectorization operator that stacks the columns of $A$ on top of one another into a vector. For two matrices $A$ and $B$, $A\otimes B$ is the Kronecker product. The set of symmetric matrices is denoted by $\mathcal{S}^{n}$. A symmetric matrix $S\in\mathcal{S}^{n}$ is said to be positive semi-definite (resp. positive definite) and written as $S\succeq 0$ (resp. $S\succ 0$) if all its eigenvalues are non-negative (resp. positive). The expression $A\succeq B$ (resp. $A\succ B$) for $A$, $B\in\mathcal{S}^{n}$ means $A-B\succeq 0$ (resp. $A-B\succ 0$). We use $\mathcal{S}_{+}^{n}$ and $\mathcal{S}_{++}^{n}$ to denote the sets of positive semi-definite matrices and positive definite matrices, respectively. For a positive semi-definite matrix $S$, $S^{\frac{1}{2}}$ is the unique positive semi-definite matrix such that $S=S^{\frac{1}{2}}S^{\frac{1}{2}}$. The $\ell_{p}$ norm for $V\in{\mathbb{R}}^{d}$ is denoted by $\|V\|_{p}=\big{(}\sum^{d}_{k=1}|V_{k}|^{p}\big{)}^{1/p}$, $p\geq 1$. The Frobenius norm and the off-diagonal $\ell_{1}$ (semi)norm of a matrix $A\in{\mathbb{R}}^{m\times n}$ are denoted by $\|A\|_{F}=\sqrt{\sum_{k=1}^{m}\sum_{l=1}^{n}|A_{kl}|^{2}}$ and $\|A\|_{1,\mathrm{off}}=\sum_{k\neq l}|A_{kl}|$, respectively. The spectral norm, i.e., the maximum singular value of a matrix $A$ is written as $\|A\|_{s}$. For a symmetric matrix A, we use $\lambda_{\max}(A)$ (resp. $\lambda_{\min}(A)$) to denote its largest (resp. smallest) eigenvalue. The maximum absolute value of all entries of a matrix $A$ is denoted by $\|A\|_{\max}$. We use $\mathcal{S}_{\mathrm{off}}^{n}$ to denote the set of $n\times n$ symmetric matrices whose diagonal elements are $0$. By $\{Y_{t,\mathrm{off}}\}_{t=1}^{T}$ we denote a sequence of symmetric matrices whose diagonal elements are $0$, i.e., $Y_{t,\mathrm{off}}\in\mathcal{S}_{\mathrm{off}}^{n}$. For a sequence of symmetric matrices $\{\Theta_{t}\}_{t=1}^{T}$, $\Theta_{uv,t}$ refers to the $(u,v)$-th element of $\Theta_{t}$, and $\Theta_{t,\mathrm{off}}$ is the copy of $\Theta_{t}$ with diagonal elements set to $0$; in particular, $\|\Theta_{t}\|_{1,\mathrm{off}}=\|\Theta_{t,\mathrm{off}}\|_{1}$. Given $\epsilon\geq 0$, we use $\mathrm{Proj}_{\cdot\succeq\epsilon I_{p}}$ to denote the projection onto $\{S\,:\,S\succeq\epsilon I_{p}\}$. The proximal operator of a function $f$ at $x$ is defined as $\mathrm{prox}_{f}(x)=\operatorname*{\arg\min}_{y}\{f(y)+\frac{1}{2}\|y-x\|^{2}\}$; for more details about the proximal operator, we refer the interested readers to Sections 12.4, 14.1, 14.2 in [1].

For a sequence of a $p$-dimensional random vector $(X_{t})$ observed at $t=1,\ldots,T$, we consider the estimation of the underlying network through the corresponding precision matrix. The latter is assumed to evolve over time and the task is to recover the break dates. More formally, we denote by $\{{\mathcal{B}}_{j}\}_{1\leq j\leq m}$ a disjoint partitioning of the set $\{1,\ldots,T\}$ such that ${\mathcal{B}}_{j}\cap{\mathcal{B}}_{j^{\prime}}=\emptyset,j\neq j^{\prime}$, $\cup_{j}{\mathcal{B}}_{j}=\{1,\ldots,T\}$ and ${\mathcal{B}}_{j}=\{T_{j-1},T_{j-1}+1,\ldots,T_{j}-1\}$. The partition of the break dates is denoted by ${\mathcal{T}}_{m}=\{T_{1}<T_{2}<\ldots<T_{m}\}$ with the convention $T_{0}=1,T_{m+1}=T+1$. Then, we assume ${\mathbb{E}}[X_{t}]=0$ and Var$(X_{t})=\Sigma_{j}$ for $t\in{\mathcal{B}}_{j}$, such that the observations indexed by elements in ${\mathcal{B}}_{j}$ are $p$-dimensional realizations of a centered random variable with variance-covariance $\Sigma_{j}$. We denote by $\Omega_{j}=\Sigma^{-1}_{j}$ the precision matrix with entries $\Omega_{uv,j},1\leq u,v\leq p$. In practice, we consider the sequence of precision matrices $\{\Theta_{1},\ldots,\Theta_{T}\}$ such that the total number of distinct matrices in the set is $m+1$ and $\Theta_{t}=\Omega_{j},\;t\in{\mathcal{B}}_{j},\;j=1,\ldots,m+1$. We are interested in estimating the unknown true number $m^{\ast}$ of unknown break dates, the true partition ${\mathcal{T}}^{\ast}_{m^{\ast}}=\{T^{\ast}_{1}<T^{\ast}_{2}<\ldots<T^{\ast}_{m^{\ast}}\}$ and the true unknown precision matrices $\Omega^{\ast}_{j}$. As a consequence, the true data generating process is assumed to be

$${\mathbb{E}}[X_{t}]=0,\;\text{Var}(X_{t})=\Sigma^{\ast}_{j},\;\Theta^{\ast}_{t}=\Omega^{\ast}_{j}=\Sigma^{\ast-1}_{j}\;\text{when}\;t=T^{\ast}_{j-1},T^{\ast}_{j-1}+1,\ldots,T^{\ast}_{j}-1,$$

and $1\leq j\leq m^{\ast}+1$, $T^{\ast}_{0}=1,T^{\ast}_{m^{\ast}+1}=T+1$ with blocks ${\mathcal{B}}^{\ast}_{j}=\{T^{\ast}_{j-1},\ldots,T^{\ast}_{j}-1\}$. While $m^{\ast}$ and the break dates are unknown, $m^{\ast}$ is typically much smaller than $T$ and, assuming the underlying network may exhibit some sparse structures, we consider the sparse estimation of $\Theta_{t}$’s and the estimation of ${\mathcal{T}}_{m}$ via a mixture of LASSO and Group Fused LASSO, which we will hereafter refer as the Group Fused D-trace LASSO (GFDtL), defined as

$$\{\widehat{\Theta}_{t}\}^{T}_{t=1}=\underset{\Theta_{t}\succeq 0,1\leq t\leq T}{\arg\;\min}\Big{\{}{\mathbb{L}}(\{\Theta_{t}\}^{T}_{t=1},\mathcal{X}_{T})+\lambda_{1}\overset{T}{\underset{t=1}{\sum}}\|\Theta_{t}\|_{1,\text{off}}+\lambda_{2}\overset{T}{\underset{t=2}{\sum}}\|\Theta_{t}-\Theta_{t-1}\|_{F}\Big{\}},$$

(2.1)

where $\mathcal{X}_{T}=(X_{1},\ldots,X_{T})$ is the sample, $\lambda_{1},\lambda_{2}$ are the tuning parameters, $\|\Theta_{t}\|_{1,\text{off}}=\sum_{k\neq l}|\Theta_{kl,t}|$ and the D-trace loss of [44] is defined as

$${\mathbb{L}}(\{\Theta_{t}\}^{T}_{t=1},\mathcal{X}_{T})=\frac{1}{T}\overset{T}{\underset{t=1}{\sum}}\Big{[}\text{tr}(\frac{1}{2}\Theta^{2}_{t}X_{t}X^{\top}_{t})-\text{tr}(\Theta_{t})\Big{]}.$$

Sparsity within the estimated precision matrix for a given block is controlled by $\lambda_{1}$, whereas $\lambda_{2}$ affects the smoothing and guarantees that the solution is piece-wise constant.

Before stating the large sample results, we define some notations and present the assumptions used hereafter. Define $\mathcal{I}^{\ast}_{j}=T^{\ast}_{j}-T^{\ast}_{j-1}$ and

$$\mathcal{I}_{\min}=\underset{1\leq j\leq m^{\ast}+1}{\min}|\mathcal{I}^{\ast}_{j}|,\;\eta_{\min}=\underset{1\leq j\leq m^{\ast}}{\min}\|\Omega^{\ast}_{j+1}-\Omega^{\ast}_{j}\|_{F},\;s^{\ast}_{\max}=\underset{1\leq j\leq m^{\ast}}{\max}\|\Omega^{\ast}_{j}\|_{F}.$$

Assumption 1-(i) relates to the properties of $(X_{t})$. Assumption 1-(ii) is a tail condition and will allow us to apply exponential inequalities for dependent processes. Assumption 2 ensures the identification of the model: it is similar to Assumption A.1 of [20] or Assumption A.2 of [30]. Assumption 3 provides conditions on $\delta_{T}$, $m^{\ast}$, $\mathcal{I}_{\min}$, $\eta_{\min}$ and the tuning parameters $\lambda_{1},\lambda_{2}$. Condition (i) concerns the convergence rate of $\delta_{T}$ to $0$. In condition (ii), the sample size in each regime may diverge with rate $T\delta_{T}$, but at a slower rate than $T$, and the number of breaks $m^{\ast}$ may diverge slowly: this is similar to Assumption A.3-(i) of [30] or Assumption H3 of [4]. It also sets the slowest rate at which $\delta_{T}$ may shrink to zero: $\delta_{T}=o(\mathcal{I}_{\min}/T)$. Conditions (iii) and (iv) specify the fastest rate at which $\delta_{T}$ may shrink to zero, which is $\delta_{T}\gg\max(\lambda_{2}/\eta_{\min},p^{2}(s^{\ast}_{\max})^{2}\log(pT)/(T\eta_{\min}^{2}))$. It is worth emphasizing that conditions (ii)-(iii) imply $p^{2}(s^{\ast}_{\max})^{2}\log(pT)=o(\mathcal{I}_{\min}\eta^{2}_{\min})$ and conditions (ii) and (iv) imply that $\lambda_{2}T=o(\eta_{\min}\mathcal{I}_{\min})$. Finally, the effect of the LASSO shrinkage through $\lambda_{1}$ does not relate to $\delta_{T}$: this is because the Group Fused LASSO penalty only allows to detect change-points. The consistency of $\widehat{T}_{j},\widehat{\Omega}_{j}$, given $\widehat{m}=m^{\ast}$, is provided in the next Theorem.

Remark Result (i) implies $\max_{1\leq j\leq m^{\ast}}T^{-1}|\widehat{T}_{j}-T^{\ast}_{j}|=o_{p}(\delta_{T})$. Since $\delta_{T}=o(1)$, this means $T^{-1}|\widehat{T}_{j}-T^{\ast}_{j}|=o_{p}(1)$. Here, $\delta_{T}$ is a key quantity to control for the rate at which $\widehat{T}_{j}/T$ converges to $T^{\ast}_{j}/T$. Note that $\delta_{T}\gg\max(\frac{\lambda_{2}}{\eta_{\min}},p^{2}(s^{\ast}_{\max})^{2}\log(pT)/(T\eta_{\min}^{2}))$, implies that the fastest convergence rate for the break ratio estimator depends on the regularization parameter $\lambda_{2}$ and $p^{2}(s^{\ast}_{\max})^{2}\log(pT)/(T\eta_{\min}^{2})$. Result (ii) relates to the consistency of the precision matrix in each regime.

It is worth noting that the true number of breaks $m^{\ast}$ is unknown. Following the common practice in the change-point literature, we assume that $m^{\ast}$ is bounded by a known conservative upper bound $m_{\max}$. Next, we define $h(A,B)\coloneqq\sup_{b\in B}\inf_{a\in A}|a-b|$ for any two sets $A$ and $B$. The next result establishes that all true break points in ${\mathcal{T}}^{\ast}_{m^{\ast}}$ can be consistently estimated by some points in $\widehat{{\mathcal{T}}}_{\widehat{m}}$.

The proof of Theorem 3.2 is done by contradiction and follows similar arguments as in the proof of Theorem 3.1. It relies on the optimality conditions from Lemma A.3. Theorem 3.2 ensures that even if the number of blocks is overestimated, there will be an estimated change-point close to each unknown true change-point.

In this section, we move on to the optimization aspects of criterion (2.1), including the dual problem, the existence of solutions, and the algorithm. Specifically, given $\mathcal{X}_{T}$, $\epsilon>0$, $\lambda_{1}>0$ and $\lambda_{2}>0$, we consider the following scaled optimization problem

$$\min_{\begin{subarray}{c}\Theta_{t}\succeq\epsilon I_{p},\\[2.84544pt] 1\leq t\leq T\end{subarray}}\left\{\sum_{t=1}^{T}\!\left[\text{tr}(\frac{1}{2}\Theta_{t}^{2}X_{t}X_{t}^{\top})\!-\!\text{tr}(\Theta_{t})\right]\!+\!\lambda_{1}T\sum_{t=1}^{T}\!\|\Theta_{t}\|_{1,\text{off}}\!+\!\lambda_{2}T\sum_{t=1}^{T-1}\!\|\Theta_{t+1}-\Theta_{t}\|_{F}\right\},$$

(4.1)

where we scaled Problem (2.1) by a factor of $T$ for numerical stability. One can also notice that in (4.1) we use $\Theta_{t}\succeq\epsilon I_{p}$ rather than $\Theta\succ 0$ as in (2.1). This choice is made for practical reasons, as setting $\epsilon>0$ ensures the non-singular solutions, and the set $\{S\,:\,S\succeq\epsilon I_{p}\}$ is closed and convex and hence the projection onto it is well defined.