Spectral Change Point Estimation for High Dimensional Time Series by Sparse Tensor Decomposition

We generally use bold-face upright letters (e.g., $\mathbf{X}$) for matrices, bold-face italic letters (e.g., $\bm{\gamma}$) for vectors, calligraphic-font letters (e.g., $\mathcal{F}$) for tensors, and superscript ^∗ to signify tensor-induced vectors (e.g., $\mathcal{F}^{*}$). For any real number $r$, $\lceil r\rceil$ denotes the ceiling function, i.e., the smallest integer greater or equal to $r$, and $\lfloor r\rfloor=-\lceil-r\rceil$ denotes the floor function. For any vector $\bm{c}\in\mathbb{R}^{p}$, $\|\bm{c}\|_{0}$ denotes its $l_{0}$ norm, and $\|\bm{c}\|$ its Euclidean norm. For any $p\times p$ real matrix $\mathbf{A}$, $\|\mathbf{A}\|$ denotes its spectral norm, and $\mathbf{A}_{i,}$ its $i$-th row. We follow the tensor convention in Kolda and Bader (2009). Specifically, for a third-order tensor $\mathcal{G}\in\mathbb{C}^{p_{1}\times p_{2}\times p_{3}}$ and some $1\leq j\leq p_{2},1\leq l\leq p_{3}$, the mode 1 fibre is given by $\mathcal{G}_{:jl}$, and the slice along mode 3 is given as $\mathcal{G}_{::l}$. When $\mathcal{G}$ is real, define the vector product with $\bm{c}^{(k)}\in\mathbb{R}^{p_{k}}$ with $k=1,2,3$, e.g. $\mathcal{G}\times_{1}\bm{c}^{(1)}\equiv\sum_{i=1}^{p_{1}}\bm{c}_{i}^{(1)}\mathcal{G}_{i::}$. For any set $\mathbb{U}$, let $|\mathbb{U}|$ be its cardinality.

The first step is to estimate the time-dependent spectral density matrix $\mathbf{f}(t,\omega)$, which requires some trade-off between spectral estimation accuracy and temporal resolution of the change point location. Similar to Adak (1998) and with no loss of generality, we partition the index set $\{1,\ldots,N\}$ into $B$ time blocks of equal length $L$, which are defined by the boundaries $\{n_{b}=Lb,b=0,...,B\}$ with $B=\lfloor N/L\rfloor$. The change point search is then on the block time scale, with the change block index set defined as $\mathbb{U}=\{u_{1\leq q\leq Q}|u_{q}=\lceil v_{q}B\rceil\}$. Note that the search can be fine-tuned using overlapping windows around any change point detected on the block scale, although we will not delve into this aspect below.

Denote $\Omega$ as a set of frequencies over which we estimate the spectrum. For instance, it can be the set of Fourier frequencies $2\pi l/L$ for $1\leq l\leq\lceil L/2\rceil$, or a subset thereof, or certain frequency band of interest, as in neuroscience. Consider the $b$-th block indexed by time $n_{b-1}+1,\ldots,n_{b}$. For any frequency $\omega\in\Omega$, define $\mathcal{F}(\omega)$ as a $p\times p\times B$ (population) spectrum tensor with $\mathcal{F}_{::b}(\omega)=\mathbf{f}_{b}(\omega)$ defined as a weighted average of $\mathbf{f}(v_{q},\omega)$ with weights proportional to number of observations within $b$-th block falling into the $q$-th segment; see the formula in (LABEL:eq:f). Then, $\hat{\mathcal{F}}(\omega)$ is defined as a $p\times p\times B$ smoothed periodogram tensor with

where $K(\cdot)$ is a kernel function (Bartlett kernel in our implementation) with bandwidth $1\leq R<L$, and $\hat{\mathbf{\Sigma}}_{b}(m)$ is the lag-$m$ autocovariance matrix for the $b$-th block, i.e.,

For a generic $p\times p\times B$ spectrum or periodogram tensor $\mathcal{G}$, define the co-spectrum CUSUM transformation $\mathcal{C}_{s,e}:\mathbb{C}^{p\times p\times B}\to\mathbb{R}^{p\times p\times(e-s)}$ on the interval $\{(s,e):1\leq s\leq e\leq B\}$ such that the $(b-s+1)$-th slice of $\mathcal{C}_{s,e}(\mathcal{G})$ for $b=s,\ldots,e-1$ is

Specifically, when $s=1$ and $e=B$, we have

Apply the $\mathcal{C}_{s,e}$ transformation to $\hat{\mathcal{F}}(\omega)$ and $\mathcal{F}(\omega)$, and for simplicity denote $\hat{\mathcal{T}}_{s,e}(\omega)=\mathcal{C}_{s,e}(\hat{\mathcal{F}}(\omega))$ and $\mathcal{T}_{s,e}(\omega)=\mathcal{C}_{s,e}(\mathcal{F}(\omega))$, and their $b$-th slices as $\hat{\mathcal{T}}_{s,b,e}(\omega)$ and $\mathcal{T}_{s,b,e}(\omega)$.

In the Supplementary Section LABEL:sec:spectheory , we provide non-asymptotic bounds on the deviation between $\hat{\mathcal{T}}_{s,e}(\omega)$ and $\mathcal{T}_{s,e}(\omega)$ uniformly over $s,e$ and $\omega$, which facilitates deriving the theoretical properties of the change point detection method.

To extract and aggregate the information about the change points as contained in the third-order tensor $\hat{\mathcal{T}}_{s,e}(\omega)$ for fixed $\omega$, we develop a projection approach to compress it into a projection vector along which the projected univariate series preserves the greatest change point signal. Consider the simple case that the spectral density function changes at some specific frequency $\omega$, across a single change point (i.e. $Q=1$). Denote $v(\omega)=v_{1}(\omega)$ as the change point location, and $u(\omega)$ the corresponding change point block index. Let $\mathbf{g}(\omega)=\mathbf{g}_{1}(\omega)=\mathrm{Re}[\mathbf{f}(v_{2},\omega)-\mathbf{f}(v_{1},\omega)]$ be the co-spectral increment.

For notation simplicity, assume that any change point is the first time point of some block (when change point is in the interior of some time blocks, just replace $u(\omega)$ below by $Bv(\omega)$), the co-spectrum CUSUM is given by

Let

and $\bm{\alpha}(\omega)=\mathrm{Norm}({\bm{\alpha}}^{\prime}(\omega))$, where $\mathrm{Norm}(\bm{c})\equiv\bm{c}/\|\bm{c}\|$ is defined as the normalization operator on vector $\bm{c}$. Let $\circ$ denotes the outer product. The eigen-decomposition of $\mathbf{g}(\omega)$ entails that

Thus, it is a special tensor in that it has a CANDECOMP/PARAFAC-decomposition whose modes 1 and 2 are sparse in the leading term, symmetric and orthogonal, and identical mode 3 across the components; see Figure 3.1 for a graphical illustration.

For the third-order tensor $\hat{\mathcal{F}}(\omega)$, we seek a projector projecting each of its frontal slice matrices into a scalar, thereby rendering the tensor into a vector for change point detection. The projected series should preserve the information about the change point as much as possible. Consider a real projector $\bm{\beta}(\omega)\in\mathbb{R}^{p}$ and $\|\bm{\beta}(\omega)\|=1$ such that the tensor $\hat{\mathcal{F}}(\omega)$ is projected into a vector $\hat{\mathcal{F}}(\omega)\times_{1}\bm{\beta}(\omega)\times_{2}\bm{\beta}(\omega)$, on which we could apply the operator $\mathcal{C}_{s,e}$ since vectors are $1\times 1\times B$ tensors. Meanwhile ${\mathcal{F}}(\omega)$ is projected into a vector ${\mathcal{F}}(\omega)\times_{1}\bm{\beta}(\omega)\times_{2}\bm{\beta}(\omega)$, whose CUSUM will achieve the maximum magnitude with the projection $\bm{\beta}(\omega)=\bm{\gamma}_{1}(\omega)$, and the maximum is

hence retaining the most signal about the change point. Thus, the oracle projection direction is the leading mode-1 component in the decomposition of the tensor $\mathcal{T}_{1,B}(\omega)$.

Below, we work on recovering the optimal projection $\bm{\gamma}_{1}(\omega)$ from the tensor $\hat{\mathcal{T}}_{1,B}(\omega)$. A state-of-art tensor decomposition technique is the alternating least squares (ALS) procedure (Kolda and Bader, 2009). This involves solving the least squares problem one tensor mode at a time, while keeping the other modes fixed, and alternating between the tensor modes. Another approach is the tensor power method, which is basically a greedy or rank-1 ALS update. This is a natural generalization of the matrix power method for extracting the dominant eigenvector, with a truncated version given by Yuan and Zhang (2013) for the sparse setting. In sparse tensor decomposition, the tensor power method can be combined with a soft thresholding operator (Allen, 2012) or a truncation step (Sun et al., 2017), and the latter provided theoretical guarantees. Since the tensor decomposition for ${\mathcal{T}}_{1,B}(\omega)$ enjoys the special features including partial orthogonality, symmetry and identical third mode, we utilize these features and develop a new Algorithm 1. It combines the tensor power method with the truncated matrix power method. For any $\bm{c}\in\mathbb{R}^{p}$, $\mathrm{Trun}(\bm{c},k)$ equals the vector resulted from retaining its largest $k$ elements in magnitude but replacing other elements of $\bm{c}$ by zero. In Algorithm 1, we first discard the boundary frontal slices of the tensor $\hat{\mathcal{T}}_{s,e}(\omega)\in\mathbb{R}^{p\times p\times(e-s)}$, for technical reasons. Specifically, let $\hat{\mathcal{T}}_{s,e}^{\dagger}(\omega)\in\mathbb{R}^{p\times p\times(e-s-2\nu)}$ be its trimmed version by discarding $\nu\asymp{R\log(Np)}/{L}\leq\log(Np)$ slices at each boundary, based on which change point detection is proceeded. For the single-change-point setting, Algorithm 1 is implemented with $s=1$ and $e=B$.

Algorithm 1 is a tensor power method, with $\bm{\alpha}^{(j)}(\omega)$ and $\bm{\gamma}^{(j)}(\omega)$ alternately updated according to lines 3 and 9. Specifically, $\bm{\gamma}^{(j)}(\omega)$ is updated by a truncated matrix power method (Yuan and Zhang, 2013), under the special case here that the matrix $\mathbf{D}^{(j)}(\omega)$ is dynamic and updated along the outer iteration. The iteration may be terminated by keeping a pre-specified maximal number of iterations, or if the changes in the $\bm{\alpha}^{(j)}(\omega)$’s and $\bm{\gamma}^{(j)}(\omega)$’s are below some threshold. For initialization $\bm{\gamma}^{(0)}(\omega)$, refer to Section LABEL:sec:initial for two approaches we provide, namely, the DSPCA approach and the truncated PCA approach. The former has theoretical guarantee and the latter is computationally efficient. Although there are two levels of iterations, the outer iteration for $\bm{\alpha}^{(j)}(\omega)$ converges very quickly thanks to the identical third mode, thus a few outer iterations would work, which we explain in Remark LABEL:rem:alpha. The convergence of $\bm{\gamma}^{(j)}(\omega)$ output by this algorithm is theoretically guaranteed by Theorem 4.1, where the benefit of exploring the sparsity structure is also illustrated.

Now, based on the CUSUM tensor $\hat{\mathcal{T}}_{s,e}(\omega)$, on interval $[s,e]$, we could obtain the estimated projector $\hat{\bm{\gamma}}_{s,e}(\omega)$, and thus project $\hat{\mathcal{T}}_{s,e}(\omega)$ into a vector. Theoretically speaking, it is preferable to implement the projector estimation and the projection-enabled change point detection with two independent co-spectrum tensors, with which the consistency of the proposed change-point detection method (see Theorem 4.2) can be established. Below, we assume the existence of two independent and identically distributed co-spectrum tensors $\hat{\mathcal{F}}(\omega)$ and $\hat{\mathcal{F}}^{\prime}(\omega)$; see Algorithm 2. Note that the output $\hat{\sigma}_{s,e}(\omega)>0$ since $\hat{\mathcal{F}}(\omega)$ is positive semidefinite, and it accounts for the variation of the projected spectral density.

Based on Algorithm 2, we could detect and estimate the single-change-point location at frequency $\omega$ by assessing the value of $|\hat{\mathcal{T}}^{*}_{s,e}(\omega)|$. In practice, the data may sustain multiple change points. Also, as previously discussed, an overall conclusion of change points across frequencies is desired. We combine the threshold $l_{1}$ aggregation (Cho and Fryzlewicz, 2015) for multiple frequencies and wild binary segmentation (Fryzlewicz, 2014) for multiple change points, which will be referred to as the wild sparsified binary segmentation. Specifically, we apply binary segmentation to isolate individual change points within random windows over which we aggregate the information across frequencies. After obtaining $\hat{\mathcal{T}}^{*}_{s,e}(\omega)$ and $\hat{\sigma}_{s,e}(\omega)$ for each Fourier frequency $\omega\in\Omega$ by Algorithm 2, let $\hat{\mathcal{F}}\in\mathbb{R}^{p\times p\times B\times|\Omega|}$ be a 4-order tensor with $\hat{\mathcal{F}}_{:::l}=\hat{\mathcal{F}}(\omega_{l})$ for $1\leq l\leq|\Omega|$. We define a thresholded sum of the CUSUM across all frequencies, as follows:

where $\mathbb{I}$ is the indicator function, and $\tau$ is the threshold. The threshold step reduces the influence of irrelevant noisy parts, and can consistently detect change points.

To implement the wild binary segmentation, we first randomly sample a large number of pairs, $(s_{1},e_{1}),\ldots,(s_{J},e_{J})$, uniformly from the set $\{(s,e):1\leq s<e\leq B\}$. Then, for each $1\leq j\leq J$, $\mathfrak{C}_{s_{j},e_{j}}(\hat{\mathcal{F}})$ is computed via (3.6). Denote the maximizer of $\mathfrak{C}_{s_{j},e_{j}}(\hat{\mathcal{F}})$ over $[{s}_{j}+\nu,{e}_{j}-\nu]$ by $\hat{u}^{[j]}$, and the corresponding maximum CUSUM value is $\mathfrak{C}_{s_{j},\hat{u}^{[j]},e_{j}}(\hat{\mathcal{F}})$. Suppose $\mathfrak{C}_{s_{j},\hat{u}^{[j]},e_{j}}(\hat{\mathcal{F}}),1\leq j\leq J$ is maximized at $j=\hat{j}$. If $\mathfrak{C}_{s_{j},\hat{u}^{[\hat{j}]},e_{j}}(\hat{\mathcal{F}})$ is not above 0, conclude that no change point is found, otherwise we declare $\hat{u}^{[\hat{j}]}$ a change point and repeat the above procedure recursively on each segment to the left and right of $\hat{u}^{[\hat{j}]}$ to find other change points. Thus, we carry out change-point detection recursively, with the search restricted to those initial time windows $(s_{j},e_{j}),1\leq j\leq J$ that lie in an interval under search. The algorithm is detailed in Algorithm 3, with the key steps in lines 6-9 depicted in Figure 3.2.

We discuss in this section some aspects about the implementation of Algorithm 3, including hyper-parameters choices and practical techniques, based on our limited experience gained from the simulation experiments reported below.

The choice of block length $L$ could be set based on the user’s preference. For example, certain neuroscience data may be segmented into uniform time intervals. By Assumption 4, $L$ should enable accurate per-block spectral density estimates while not containing more than one change point. In practice, we recommend setting $L$ to be between 50 and 100. Since $R\geq(L/\log(Np))^{1/3}$ by Assumption 5 (iii), we find $R=\lfloor L^{1/3}\rfloor\lor 1$ works well for our change point detection. For Algorithms 1 and 3, we specify a trimming distance $\nu\asymp{R\log(Np)}/{L}\geq(\log(Np)/{L})^{2/3}$ for technical reasons. In practice, we set $\nu=0$ in Algorithm 1 to exploit full information in finding the projector, and $\nu=\lfloor(B\log(Np))^{2/3}/15\rfloor\lor 1$ in Algorithm 3. Also, to avoid spurious change point detection, we require $\mathfrak{C}_{s_{j},e_{j}}(\hat{\mathcal{F}})$ to stay positive in a small section around the detected change point. Specifically, let $\hat{u}^{[j]}=\arg\max_{u}\mathfrak{C}_{s_{j},u,e_{j}}(\hat{\mathcal{F}})$ within the set $\{u:\min(u-s_{j},e_{j}-u)>\nu,\mathfrak{C}_{s_{j},b,e_{j}}(\hat{\mathcal{F}})>0,|b-u|<\nu/4\}$.

As for the sparsity $k$, we use a heuristic approach to determine an upper bound of the number of series components sustaining at least one change. We do this by implementing Algorithm 3 with each univariate component time series of $\mathbf{X}$, where projection is not needed. Then, we set $k$ to be the number of such time series found to have some change points. The series-by-series approach tends to find spurious change points for stationary time series, making it non-competitive as a standalone method but providing a reasonable choice of $k$ for our proposed method since Theorem 4.2 requires $k$ to be greater than $k_{0}$. We recommend the wild interval number $J=500$. In the simulation studies, we found that even setting $J=30$ yielded satisfactory results.