Matrix factorization for multivariate time series analysis

TL;DR

This paper investigates the statistical performance of matrix factorization in multivariate time series analysis, extending existing results to time-dependent data and showing improved convergence rates with structured series.

Contribution

It extends matrix estimation theory from i.i.d. data to time series and demonstrates enhanced convergence rates for structured series like periodic or smooth ones.

Findings

Extended matrix factorization results to time series data.

Proved improved convergence rates for structured time series.

Provided theoretical insights into statistical performance.

Abstract

Matrix factorization is a powerful data analysis tool. It has been used in multivariate time series analysis, leading to the decomposition of the series in a small set of latent factors. However, little is known on the statistical performances of matrix factorization for time series. In this paper, we extend the results known for matrix estimation in the i.i.d setting to time series. Moreover, we prove that when the series exhibit some additional structure like periodicity or smoothness, it is possible to improve on the classical rates of convergence.

Equations

$$\mathbf{X}=(x_{i,t}{)}_{(i,t)\in [[1,d]]\times [[1,T]]}.$$

$$\mathbf{X}=\mathbf{M}+\epsilon ,$$

$${\mathbf{W}}_{i,t}=\sum _{n=0}^{\infty }{\mathbf{U}}_{i,n}{\mathbf{e}}_n\left(\frac{t}{T}\right)$$

$${\mathbf{W}}_{i,t}\simeq \sum _{n=0}^{\tau -1}{\mathbf{U}}_{i,n}{\mathbf{e}}_n\left(\frac{t}{T}\right)$$

$$\mathbf{X}=\mathbf{M}+\epsilon =\mathbf{UV\Lambda }+\epsilon ,$$

$$\Vert \mathbf{A}{\Vert }_F=\text{trace}(\mathbf{A}{\mathbf{A}}^{\ast }{)}^{1/2}.$$

$$\Vert \mathbf{A}{\Vert }_{\text{op}}=\underset{\Vert x\Vert =1}{\sup }\Vert \mathbf{A}x\Vert$$

$$\tilde{\mathbf{X}}=\tilde{\mathbf{M}}+\tilde{\epsilon }$$

$${\widehat{\tilde{\mathbf{M}}}}_{\mathcal{S}}\in \arg \underset{\mathbf{A}\in \mathcal{S}}{\min }\tilde{r}(\mathbf{A})$$

$$\tilde{r}(\mathbf{A})=\Vert \mathbf{A}-\tilde{\mathbf{X}}{\Vert }_F^2\text{ ; }\forall \mathbf{A}\in {\mathcal{M}}_{d,\tau }(\mathbb{R}).$$

$${\Vert {\widehat{\tilde{\mathbf{M}}}}_{\mathcal{S}}-\tilde{\mathbf{M}}\Vert }_F^2.$$

$${\Vert {\hat{\mathbf{M}}}_{\mathcal{S}}-\mathbf{M}\Vert }_F^2$$

$${\mathfrak{K}}_{\epsilon }:=\underset{\Vert x\Vert =1}{\sup }\underset{p\in [1,\infty [}{\sup }{p}^{-1/2}\mathbb{E}(|⟨{\epsilon }_{1,.}{\Sigma }_{\epsilon }^{-1/2},x⟩{|}^p{)}^{1/p}<\infty .$$

$$\underset{\Vert x\Vert =1}{\sup }\underset{p\in [1,\infty [}{\sup }{p}^{-1/2}\mathbb{E}(|⟨X,x⟩{|}^p{)}^{1/p}=C$$

$$\frac{1}{dT}{\Vert {\hat{\mathbf{M}}}_{\mathcal{S}}-\mathbf{M}\Vert }_F^2⩽\frac{1+\lambda }{1-\lambda }\cdot \underset{\mathbf{A}\in \mathcal{S}}{\min }\frac{1}{dT}{\Vert \mathbf{A}-\mathbf{M}\Vert }_F^2+\frac{4\mathfrak{c}{\tilde{\mathfrak{K}}}_{\epsilon }\Vert {\Sigma }_{\epsilon }{\Vert }_{\text{op}}}{\lambda (1-\lambda )}\cdot \frac{k(d+T+s)}{dT}$$

$$\frac{1}{dT}{\Vert {\hat{\mathbf{M}}}_{\mathcal{S}}-\mathbf{M}\Vert }_F^2⩽\frac{4\mathfrak{c}{\tilde{\mathfrak{K}}}_{\epsilon }\Vert {\Sigma }_{\epsilon }{\Vert }_{\text{op}}}{\lambda (1-\lambda )}\cdot \frac{k(d+T+s)}{dT}.$$

$$\Vert {\Sigma }_{\epsilon }{\Vert }_{\text{op}}={\sigma }^2.$$

$$\Vert {\Sigma }_{\epsilon }{\Vert }_{\text{op}}=\underset{t\in [[1,T]]}{\max }\mathrm{var}({\epsilon }_{1,t}).$$

$${\Sigma }_{\epsilon }={\sigma }^2\left(\begin{array}{ccccc}1+{\theta }^2& -\theta & 0& \dots & 0\\ -\theta & 1+{\theta }^2& 0& \dots & 0\\ ⋮& & \ddots & & ⋮\\ 0& \dots & 0& 1+{\theta }^2& -\theta \\ 0& \dots & 0& -\theta & 1+{\theta }^2\end{array}\right)$$

$$\Vert {\Sigma }_{\epsilon }{\Vert }_{\text{op}}={\sigma }^2\left[1+{\theta }^2-2\theta \underset{\ell \in [[1,T]]}{\min }\cos \left(\frac{\ell \pi }{1+T}\right)\right]⩽{\sigma }^2(1+\theta {)}^2.$$

$${\Sigma }_{\epsilon }={\sigma }^2\left(\begin{array}{ccccc}1& \rho & {\rho }^2& \dots & {\rho }^{T-1}\\ \rho & 1& \rho & \dots & {\rho }^{T-2}\\ ⋮& & \ddots & & ⋮\\ {\rho }^{T-2}& \dots & \rho & 1& \rho \\ {\rho }^{T-1}& \dots & {\rho }^2& \rho & 1\end{array}\right)={\sigma }^2\left[{\mathbf{I}}_T+\sum _{t=1}^{T-1}{\rho }^t\left({\mathbf{J}}_T^t+({\mathbf{J}}_T^{\ast }{)}^t\right)\right].$$

$${\mathbf{J}}_T=\left(\begin{array}{ccccc}0& 1& 0& \dots & 0\\ 0& 0& 1& \dots & 0\\ ⋮& & \ddots & \ddots & ⋮\\ 0& \dots & 0& 0& 1\\ 0& \dots & 0& 0& 0\end{array}\right).$$

$$\Vert {\Sigma }_{\epsilon }{\Vert }_{\mathrm{op}}⩽{\sigma }^2\left(1+2\sum _{t=1}^T|\rho {|}^t\right)⩽{\sigma }^2\left(1+\frac{2|\rho |}{1-|\rho |}\right)={\sigma }^2\frac{1+|\rho |}{1-|\rho |}.$$

$${\Vert {\widehat{\tilde{\mathbf{M}}}}_{\mathcal{S}}-\tilde{\mathbf{M}}\Vert }_F^2⩽\frac{1+\lambda }{1-\lambda }\cdot \underset{\mathbf{A}\in \mathcal{S}}{\min }{\Vert \mathbf{A\Lambda }-\tilde{\mathbf{M}}\Vert }_F^2+\frac{4\mathfrak{c}k}{\lambda (1-\lambda )}(d+\tau +s){\tilde{\mathfrak{K}}}_{\tilde{\epsilon }}\Vert {\Sigma }_{\epsilon {\mathbf{\Lambda }}^{+}}{\Vert }_{\text{op}}$$

$$\Vert {\Sigma }_{{\epsilon }_{1,.}\mathbf{C}}{\Vert }_{\text{op}}⩽\Vert {\Sigma }_{\epsilon }{\Vert }_{\text{op}}\Vert {\mathbf{C}}^{\ast }\mathbf{C}{\Vert }_{\text{op}}.$$

$${\Vert {\hat{\mathbf{M}}}_{\mathcal{S}}-\mathbf{M}\Vert }_F^2={\Vert ({\widehat{\tilde{\mathbf{M}}}}_{\mathcal{S}}-\tilde{\mathbf{M}})\mathbf{\Lambda }\Vert }_F^2⩽{\Vert {\widehat{\tilde{\mathbf{M}}}}_{\mathcal{S}}-\tilde{\mathbf{M}}\Vert }_F^2\Vert \mathbf{\Lambda }{\mathbf{\Lambda }}^{\ast }{\Vert }_{\mathrm{op}}$$

$${\Vert \mathbf{A}-\tilde{\mathbf{M}}\Vert }_F^2={\Vert (\mathbf{A\Lambda }-\mathbf{M}){\mathbf{\Lambda }}^{+}\Vert }_F^2⩽{\Vert \mathbf{A\Lambda }-\mathbf{M}\Vert }_F^2\Vert \mathbf{\Lambda }{\mathbf{\Lambda }}^{+}\Vert (\mathbf{\Lambda }{\mathbf{\Lambda }}^{\ast }{)}^{-1}{\Vert }_{\mathrm{op}}.$$

$$\frac{1}{dT}{\Vert {\hat{\mathbf{M}}}_{\mathcal{S}}-\mathbf{M}\Vert }_F^2⩽\frac{1+\lambda }{1-\lambda }\cdot \underset{\mathbf{A}\in \mathcal{S}}{\min }\frac{1}{dT}{\Vert \mathbf{A\Lambda }-\mathbf{M}\Vert }_F^2+\frac{4\mathfrak{c}{\tilde{\mathfrak{K}}}_{\epsilon }\Vert {\Sigma }_{\epsilon }{\Vert }_{\text{op}}}{\lambda (1-\lambda )}\cdot \frac{k(d+\tau +s)}{dT}$$

$$\mathbf{\Lambda }=({\mathbf{I}}_{\tau }|\dots |{\mathbf{I}}_{\tau })\in {\mathcal{M}}_{\tau ,T}(\mathbb{R}).$$

$$\mathbf{\Lambda }{\mathbf{\Lambda }}^{\ast }=\frac{T}{\tau }{\mathbf{I}}_{\tau }\Rightarrow \Vert \mathbf{\Lambda }{\mathbf{\Lambda }}^{\ast }{\Vert }_{\mathrm{op}}=\frac{T}{\tau }\text{ and }\Vert (\mathbf{\Lambda }{\mathbf{\Lambda }}^{\ast }{)}^{-1}{\Vert }_{\text{op}}=\frac{\tau }{T}.$$

Full text

Acknowledgements. The authors gratefully acknowledge Maxime Ossonce for discussions which helped to improve the paper. The first author was working at CREST, ENSAE Paris when this paper was written; he gratefully acknowledges financial support from Labex ECODEC (ANR-11-LABEX-0047).

1. Introduction

Matrix factorization is a very powerful tool in statistics and data analysis. It was used as early as in the 70’s in econometrics in reduced-rank regression [27, 22, 29]. There, matrix factorization is mainly a tool to estimate a coefficient matrix under a low-rank constraint. There was recently a renewed interest in matrix factorization as a data analysis tool for huge datasets. Nonnegative matrix factorization (NMF) was introduced by [35] as a tool to represent a huge number of objects as linear combinations of elements of “parts” of objects. The method was indeed applied to large facial image datasets and the dictionary indeed contained typical parts of faces. Since then, various methods of matrix factorization were successfully applied such various fields as collaborative filtering and recommender systems on the Web [34, 51], document clustering [46], separation of sources in audio processing [42], missing data imputation [26], quantum tomography [23, 25, 53, 7, 39], medical image processing [21] topics extraction in texts [43] or transports data analysis [12]. Very often, matrix factorization provides interpretable and accurate representations of the data matrix as the product of two much smaller matrices. The theoretical performances of matrix completion were studied in a series of papers by Candès with many co-authors [10, 11, 9]. Minimax rates for matrix completion and more general matrix estimation problems were derived in [32, 8, 30, 31, 41]. Bayesian estimators and aggregation procedures were studied in [1, 47, 38, 2, 37, 4, 17, 36, 16, 15].

To apply matrix factorization techniques to multivariate time series is a very natural idea. First, the low-rank structure induced by the factorization leads to high correlations that are indeed observed in some applications (this structure is actually at the core of cointegration models in econometrics [20, 28, 5]). Moreover, the factorization provides a decomposition of each series in a dictionary which member that can be interpreted as latent factors used for example in state-space models, see e.g. Chapter 3 in [33]. For this reasons, matrix factorization was used in multivariate time series analysis beyond econometrics: electricity consumptions forecasting [18, 40], failure detection in transports systems [48], collaborative filtering [24], social media analysis [45] to name a few.

It is likely that the temporal structure in the data can be exploited to obtain an accurate and sensible factorization: autocorrelation, smoothness, periodicity… Indeed, while some authors use matrix factorization as a black box for data analysis, others propose in a way or another to adapt the algorithm to the temporal structure of the data [54, 45, 14, 24]. However, there is no theoretical guarantee that this leads to better predictions or better rates of convergence. Moreover, the aforementionned theoretical studies [10, 11, 9, 32, 8, 30, 31, 41] all assumed i.i.d noise, strongly limiting their applicability to study algorithms designed for time series such as in [54]. The objective of this paper is to address both issues.

Consider for example that one observes a $d$ series $(x_{i,t})_{t=1}^{T}=\mathbf{X}$ and assume that $\mathbf{X}=\mathbf{M}+\varepsilon$ where $\mathbf{M}$ is a rank $k$ matrix and $\varepsilon$ is some noise. In a first time, assume that entries $\varepsilon_{i,t}$ of $\varepsilon$ are i.i.d with variance $\sigma^{2}$. Theorem 3 in [32] implies that there is an estimator $\hat{\mathbf{M}}_{1}$ of $\mathbf{M}$, such that $\frac{1}{dT}\|\hat{\mathbf{M}}_{1}-\mathbf{M}\|_{F}^{2}=\mathcal{O}(\sigma^{2}\frac{k(d+T)}{dT})$, up to log terms. Moreover, Theorem 5 in the same paper shows that this rate cannot be improved. Here, we propose an estimator $\hat{\mathbf{M}}=\hat{\mathbf{U}}\hat{\mathbf{V}}$, where $\hat{\mathbf{U}}$ is a $d\times k$ matrix and $\hat{\mathbf{V}}$ is $k\times T$. We study this estimator under the assumption that the rows $\varepsilon_{i,\cdot}$ of $\varepsilon$ are independent, centered, with covariance matrix $\Sigma_{\varepsilon}$, allowing a temporal dependence in the noise. We prove that $\frac{1}{dT}\|\hat{\mathbf{M}}-\mathbf{M}\|_{F}^{2}=\mathcal{O}(\left\|\Sigma_{\varepsilon}\right\|_{{\rm op}}\frac{k(d+T)}{dT})$ where $\|\Sigma_{\varepsilon}\|_{{\rm op}}$ is the operator norm of $\Sigma_{\varepsilon}$. Note that in the i.i.d case $\Sigma_{\varepsilon}=\sigma^{2}\mathbb{I}_{T}$, we recover the rate of [32] as $\|\Sigma_{\varepsilon}\|_{{\rm op}}=\sigma^{2}$. However, our result is more general: we provide examples where the noise is non i.i.d and we still have a control on $\|\Sigma_{\varepsilon}\|_{{\rm op}}$. For example, when the noise is row-wise AR(1), that is $\varepsilon_{i,t+1}=\rho\varepsilon_{i,t}+\eta_{i,t}$ where the $\eta_{i,t}$ are i.i.d with variance $\sigma^{2}$ and $|\rho|<1$, we have $\|\Sigma_{\varepsilon}\|_{{\rm op}}\leqslant\sigma^{2}\frac{1+|\rho|}{1-|\rho|}$. Moreover, our estimator can be tuned to take into account a possible periodicity or smoothness of the series. This is done by rewriting $\hat{\mathbf{W}}=\hat{\mathbf{V}}\mathbf{\Lambda}$ where $\mathbf{\Lambda}$ is a $\tau\times T$ matrix encoding the temporal structure, and $\tau\leqslant T$. In this case, we always improve on the rate $\mathcal{O}(\left\|\Sigma_{\varepsilon}\right\|_{{\rm op}}\frac{k(d+T)}{dT})$.

We obtain the following rates, for some constant $C(\beta,L)$:

[TABLE]

All the results are first stated under a known structure, that is, we assume that we know the rank $k$, the period $\tau$ or the smoothness $\beta$ of the series. We provide at the end of the paper a model selection procedure that allows to obtain the same rates of convergence without assuming this prior knowledge.

Finally, we should mention the nice paper [44] where the authors studied time-evolving adjacency matrices for graphs with autoregressive features. However, the rows of an adjacency matrix are not interpreted as time series, so the objective of this work is quite different from ours.

The paper is organized as follows. In Section 2 we introduce the notations that will be used in all the paper. In Subection 3.1 we study matrix factorization without additional temporal structure. In Section 3.2, we study the estimator $\hat{\mathbf{M}}=\hat{\mathbf{U}}\hat{\mathbf{V}}\mathbf{\Lambda}$ in the general case, and show how it improves the rates of convergence for a well chosen matrix $\mathbf{\Lambda}$ for periodic and/or smooth series. Finally, adaptation to unknown rank, periodicity and/or smoothness is tackled in Section 4. The proofs are given in Section 5.

2. Setting of the problem and notation

Assume that we observe a multivariate series

[TABLE]

where $d\in\mathbb{N}^{*}$ and $T\in\mathbb{N}\backslash\{0,1\}$. This multivariate series is modelled as a stochastic process. We actually assume that

[TABLE]

where $\varepsilon$ is a noise and $\mathbf{M}$ is a matrix of rank $k\in\llbracket 1,T\rrbracket$. Then, there exist $\mathbf{U}\in\mathcal{M}_{d,k}(\mathbb{R})$ and $\mathbf{W}\in\mathcal{M}_{k,T}(\mathbb{R})$ such that $\mathbf{M}=\mathbf{U}\mathbf{W}$. We will refer to $\mathbf{W}$ as the dictionary or as the latent series.

We also want to model more structure in $\mathbf{M}$. This is done by rewritting $\mathbf{W}=\mathbf{V}\mathbf{\Lambda}$, where $\tau\in\mathbb{N}\setminus\{0\}$, $\mathbf{V}\in\mathcal{M}_{k,\tau}$ and $\mathbf{\Lambda}\in\mathcal{M}_{\tau,T}$, where $\mathbf{\Lambda}$ is a known matrix. The matrix $\Lambda$ depends on the structure assumed on $\mathbf{M}$.

Example 2.1** (Periodic series).**

Assume that $T=p\tau$ with $p\in\mathbb{N}^{*}$ for the sake of simplicity. To assume that the latent series in the dictionary $\mathbf{W}$ are $\tau$-periodic is exactly equivalent to writing $\mathbf{W}=\mathbf{V}\mathbf{\Lambda}$ where $\mathbf{V}\in\mathcal{M}_{k,\tau}(\mathbb{R})$ and $\mathbf{\Lambda}=(\mathbf{I}_{\tau}|\dots|\mathbf{I}_{\tau})\in\mathcal{M}_{\tau,T}(\mathbb{R})$ is defined by blocks, $\mathbf{I}_{\tau}$ being the indentity matrix in $\mathcal{M}_{\tau,\tau}(\mathbb{R})$.

Example 2.2** (Smooth series).**

We can assume that the series in $\mathbf{W}$ are smooth. For example, say that they belong a a Sobolev space with smoothness $\beta$, we have

[TABLE]

where $(\mathbf{e}_{n})_{n\in\mathbb{N}}$ is the Fourier basis (the definition of a Sobolev space is reminded is Section 3.2 below). Of course, there are infinitely many coefficients $\mathbf{U}_{i,n}$ and to estimate them all is not feasible, however, for $\tau$ large enough, the approximation

[TABLE]

will be suitable, and can be rewritten as $\mathbf{W}=\mathbf{U}\mathbf{\Lambda}$ where $\mathbf{\Lambda}_{i,t}=e_{i}(t/T)$. More details will be given in Section 3.2, where we actually cover more general basis of functions.

So our complete model will finally be written as

[TABLE]

where $\mathbf{U}\in\mathcal{M}_{d,k}(\mathbb{R})$ and $\mathbf{V}\in\mathcal{M}_{k,\tau}(\mathbb{R})$ are unknown, but $\tau\leqslant T$ and $\mathbf{\Lambda}\in\mathcal{M}_{\tau,T}(\mathbb{C})$ such that ${\rm rank}(\mathbf{\Lambda})=\tau$ are known (note that the unstructured case corresponds to $\tau=T$ and $\mathbf{\Lambda}=\mathbf{I}_{T}$).

Note that more constraint can be imposed on the estimator. For example, in nonnegative matrix factorization [35], one imposes that all the entries in $\hat{\mathbf{U}}$ and $\hat{\mathbf{W}}$ are nonnegative. Here, we will more generally assume that $\hat{\mathbf{U}}\hat{\mathbf{V}}$ belongs to some prescribed subset $\mathcal{S}\subseteq\mathcal{M}_{d,T}(\mathbb{R})$.

In what follows, we will consider two norms on $\mathcal{M}_{d,T}$. For a matrix $A$, the Frobenius norms is given by

[TABLE]

and the operator norm by

[TABLE]

where $\|\cdot\|$ is the Euclidean norm on $\mathbb{R}^{T}$.

2.1. Estimation by empirical risk minimization

By multiplying both sides in (2) by the pseudo-inverse $\mathbf{\Lambda}^{+}=\mathbf{\Lambda}^{*}(\mathbf{\Lambda}\mathbf{\Lambda}^{*})^{-1}$, we obtain the “simplified model”

[TABLE]

with $\widetilde{\mathbf{X}}=\mathbf{X}\mathbf{\Lambda}^{+}$, $\widetilde{\mathbf{M}}=\mathbf{UV}$ and $\widetilde{\varepsilon}=\varepsilon\mathbf{\Lambda}^{+}$. In this model, the estimation of $\widetilde{\mathbf{M}}$ can be done by empirical risk minimization:

[TABLE]

where

[TABLE]

Therefore, we can define the estimator $\widehat{\mathbf{M}}_{\mathcal{S}}=\widehat{\widetilde{\mathbf{M}}}_{\mathcal{S}}\mathbf{\Lambda}$ of $\mathbf{M}$.

In Section 3, we study the statistical performances of this estimator. The first step is done in Subsection 3.1, where we derive upper bounds on

[TABLE]

The corresponding upper bounds on

[TABLE]

are derived in Subsection 3.2.

3. Oracle inequalities

Throughout this section, assume that $\varepsilon$ fulfills the following..

Assumption 3.1**.**

The rows of $\varepsilon$ are independent and have the same $T$-dimensional sub-Gaussian distribution, with second moment matrix $\Sigma_{\varepsilon}$. Moreover, $\varepsilon_{1,.}\Sigma_{\varepsilon}^{-1/2}$ is isotropic and has a finite sub-Gaussian norm

[TABLE]

In the sequel, we also consider $\widetilde{\mathfrak{K}}_{\varepsilon}:=\mathfrak{K}_{\varepsilon}^{2}\vee\mathfrak{K}_{\varepsilon}^{4}$.

We remind (see e.g Chapter 1 in [13]) that when $X\sim\mathcal{N}(0,\mathbf{I}_{n})$,

[TABLE]

for some universal constant $C>0$ (that is, $C$ does not depend on $n$). Thus, for Gaussian noise, Assumption 3.1 is satisfied and $\mathfrak{K}_{\varepsilon}=C$ does not depend on the dimension $T$.

3.1. The case $\mathbf{\Lambda}=\mathbf{I}_{T}$

In this subsection only, we assume that $\mathbf{\Lambda}=\mathbf{I}_{T}$ (and thus $\tau=T$). So the simplified model is actually the original model $\widetilde{\mathbf{X}}=\mathbf{X}$, $\widetilde{\mathbf{M}}=\mathbf{M}$, $\widetilde{\mathbf{\varepsilon}}=\mathbf{\varepsilon}$ and $\widehat{\widetilde{\mathbf{M}}}_{\mathcal{S}}=\widehat{\mathbf{M}}_{\mathcal{S}}$.

Theorem 3.2**.**

Under Assumption 3.1, for every $\lambda\in]0,1[$ and $s\in\mathbb{R}_{+}$,

[TABLE]

with probability larger than $1-2e^{-s}$.

As a consequence, if we have indeed $\mathbf{M}\in\mathcal{S}$, then with large probability,

[TABLE]

Thus, we recover the rate $\mathcal{O}(\|\Sigma_{\varepsilon}\|_{{\rm op}}\frac{k(d+T)}{dT})$ claimed in the introduction.

Remark 3.3**.**

Since the bound relies on the constant $\|\Sigma_{\varepsilon}\|_{{\textrm{op}}}$, let us provide its value in some special cases:

1. (1)

If ${\rm cov}(\varepsilon_{1,t},\varepsilon_{1,t^{\prime}})=\sigma^{2}\mathbf{1}_{\{t=t^{\prime}\}}$ then

[TABLE]

More generally, when $\varepsilon_{1,1},\dots,\varepsilon_{1,T}$ are uncorrelated,

[TABLE] 2. (2)

Let $(\eta_{t})_{t\in\mathbb{Z}}$ be a white noise of standard deviation $\sigma>0$ and assume that there exists $\theta\in\mathbb{R}^{*}$ such that $\varepsilon_{1,t}=\eta_{t}-\theta\eta_{t-1}$ for every $t\in\llbracket 1,T\rrbracket$. In other words, $(\varepsilon_{1,t})_{t=1,\dots,T}$ is the restriction of a MA(1) process to $\llbracket 1,T\rrbracket$. So,

[TABLE]

and then

[TABLE] 3. (3)

Let $(\eta_{t})_{t\in\mathbb{Z}}$ be a white noise of standard deviation $\sigma>0$ and assume that there is a $\rho$ with $|\rho|<1$ such that $\varepsilon_{1,t}=\rho\varepsilon_{1,t-1}+\eta_{t}$. So $(\varepsilon_{1,t})_{t=1,\dots,T}$ is the restriction of a AR(1) process to $\llbracket 1,T\rrbracket$. So,

[TABLE]

where

[TABLE]

As $\|\mathbf{J}_{T}\|_{\rm op}=1$, we have

[TABLE]

3.2. The general case

Let us now come back to the general case. An application of Theorem 3.2 to the “simplified model” (2.1) shows that for any $\lambda\in]0,1[$ and $s\in\mathbb{R}_{+}$,

[TABLE]

with probability larger than $1-2e^{-s}$.

In order to obtain the desired bound on $\|\widehat{\mathbf{M}}_{\mathcal{S}}-\mathbf{M}\|_{F}^{2}$, we must now understand the behaviour of $\|\Sigma_{\varepsilon\mathbf{\Lambda}^{+}}\|_{{\textrm{op}}}$ and $\widetilde{\mathfrak{K}}_{\widetilde{\varepsilon}}$.

Lemma 3.4**.**

For any matrix $\mathbf{C}\in\mathcal{M}_{T,\tau}(\mathbb{C})$,

[TABLE]

The situation regarding $\widetilde{\mathfrak{K}}_{\widetilde{\varepsilon}}=\mathfrak{K}_{\widetilde{\varepsilon}}^{2}\vee\mathfrak{K}_{\widetilde{\varepsilon}}^{4}$ is different, we are not aware of a general simple upper bound on $\mathfrak{K}_{\tilde{\varepsilon}}=\mathfrak{K}_{\varepsilon\mathbf{\Lambda}^{+}}$ in terms of $\mathfrak{K}_{\varepsilon}$ and $\mathbf{\Lambda}^{+}$. Still, there are two cases where we actually have $\mathfrak{K}_{\tilde{\varepsilon}}=\mathfrak{K}_{\varepsilon}$. Indeed, in the Gaussian case, $\mathfrak{K}_{\tilde{\varepsilon}}=\mathfrak{K}_{\varepsilon}=C$, see (4) above. For non Gaussian noise, we have the following result.

Lemma 3.5**.**

Assume that there is $c(\tau,T)>0$ such that $\mathbf{\Lambda}\mathbf{\Lambda}^{*}=c(\tau,T)\mathbf{I}_{\tau}$. If $\Sigma_{\varepsilon}=\sigma^{2}\mathbf{I}_{T}$ with $\sigma>0$, then $\mathfrak{K}_{\widetilde{\varepsilon}}=\mathfrak{K}_{\varepsilon}$.

Note that the assumption on $\mathbf{\Lambda}$ is fullfilled by the examples covered in Subsections 3.3 and (3.4).

The previous discussion legitimates the following assumption.

Assumption 3.6**.**

$\mathfrak{K}_{\widetilde{\varepsilon}}\leqslant\mathfrak{K}_{\varepsilon}$.

Finally, note that

[TABLE]

and in the same way

[TABLE]

By Inequality (5) together with Lemmas 3.4 and 3.5, we obtain the following result.

Corollary 3.7**.**

Fix $\lambda\in]0,1[$ and $s\in\mathbb{R}_{+}$. Under Assumption 3.1 and Assumption 3.6,

[TABLE]

with probability larger than $1-2e^{-s}$.

Corollary 3.7 provides an oracle inequality: it says that our estimator provides the optimal tradeoff between a variance term in $\|\Sigma_{\varepsilon}\|_{{\textrm{op}}}k(d+\tau)/(dT)$, and a bias term. The bias term is the distance of $\mathbf{M}$ to its best approximation by a matrix of the form $\mathbf{A\Lambda}$. In order to explicit the rates of convergence, assumptions can be made to upper-bound the bias term. We now apply Corollary 3.7 in the case of periodic time series, and then in the case of smooth time series. In each case, we explicit the bias term and the rate of convergence.

3.3. Application: periodic time series

In the case of $\tau$-periodic time series, remind that we assumed for simplicity that there is an integer $p$ such that $\tau p=T$ and we defined

[TABLE]

Then

[TABLE]

Therefore, by Corollary 3.7, for every $\lambda\in]0,1[$ and $s\in\mathbb{R}_{+}$, under Assumptions 3.1 and 3.6,

[TABLE]

with probability larger than $1-2e^{-s}$. Now, define

[TABLE]

and assume that $\mathbf{M}\in\mathcal{S}$. Then,

[TABLE]

which is indeed an improvement with respect to the rate obtained without taking the periodicity into account, that is $\mathcal{O}(\|\Sigma_{\varepsilon}\|_{{\rm op}}\frac{k(d+T)}{dT})$.

3.4. Application: time series with smooth trend

Assume we are given a dictionary of functions $(\mathbf{e}_{n})_{|n|\leqslant N}$ for some finite $N\in\mathbb{N}$. This dictionary can for example be a finite subset of a basis of an Hilbert space $(\mathbf{e}_{n})_{n\in\mathbb{Z}}$, like the Fourier basis or a wavelet basis.

Define

[TABLE]

Note that $\mathbf{\Lambda}_{N}$ is a $\tau\times T$ matrix where $\tau=2N+1$.

Assume that

[TABLE]

This implies that $\|(\mathbf{\Lambda}_{N}\mathbf{\Lambda}_{N}^{*})^{-1}\|_{{\textrm{op}}}=1/T$ and $\|\mathbf{\Lambda}_{N}\mathbf{\Lambda}_{N}^{*}\|_{{\textrm{op}}}=T$. This can be the case for a well-chosen basis, otherwise, we can apply the Gram-Schmidt to the dictionary of functions.

Example 3.8**.**

(Fourier’s basis) Consider the Fourier basis $(\mathbf{e}_{n})_{n\in\mathbb{Z}}$ defined by

[TABLE]

On the one hand, for every $n\in\llbracket-N,N\rrbracket$ and $t\in\llbracket 1,T\rrbracket$, $|e_{n}(t/T)|=1$. On the other hand, for every $m,n\in\llbracket-N,N\rrbracket$ such that $m\not=n$,

[TABLE]

Therefore, by Corollary 3.7, for every $\lambda\in]0,1[$ and $s\in\mathbb{R}_{+}$, under Assumptions 3.1 and 3.6,

[TABLE]

with probability larger than $1-2e^{-s}$. We will now show the consequences of these results when the rows of $\mathbf{M}$ are smooth in the sense that they belong to a given Sobolev ellipsoid. In this case, we will not have a $\mathbf{A}$ such that $\left\|\mathbf{A}\mathbf{\Lambda}-\mathbf{M}\right\|_{F}^{2}=0$, but this quantity will be small and can be controlled as a function of $N$. We introduce a few definitions.

Definition 3.9**.**

The Sobolev ellipsoid $W(\beta,L)$ is the set of functions $f:[0,1]\rightarrow\mathbb{R}$ such that $f$ is $\beta-1$ times differentiable, $f^{(\beta-1)}$ is absolutely continuous and

[TABLE]

From now, we assume that $\mathbf{e}_{n}(x)=e^{2i\pi nx}$ is the Fourier basis. It is well-known from Chapter 1 in [50] that any $f\in W(\beta,L)$ and $x\in[0,1]$,

[TABLE]

and that there is a (known) constant $C(\beta,L)>0$ such that

[TABLE]

Definition 3.10**.**

We define $\mathcal{S}(k,\beta,L)\subset\mathcal{M}_{d,T}(\mathbb{R})$ as the set of matrices $\mathbf{M}$ such that $\mathbf{M}=\mathbf{U}\mathbf{W}$, $\mathbf{U}\in\mathcal{M}_{k,T}(\mathbb{R})$, $\mathbf{W}\in\mathcal{M}_{d,k}(\mathbb{R})$ and

1. (1)

For any $i\in\llbracket 1,d\rrbracket$, $\|\mathbf{U}_{i,\cdot}\|^{2}\leqslant 1$, 2. (2)

For any $\ell\in\llbracket 1,k\rrbracket$ and $t\in\llbracket 1,T\rrbracket$, $\mathbf{W}_{\ell,t}=f_{\ell}\left(\frac{t}{T}\right)$ for some $f_{\ell}\in W(\beta,L)$.

Denote $\mathbf{V}_{N,\mathbf{W}}=\left(c_{n}(f_{\ell})\right)_{\ell\leqslant k,|n|\leqslant N}$.

Then (7) implies

[TABLE]

Pluging this into (6) gives

[TABLE]

If $\beta$ is known, an adequate optimization with respect to $N$ gives the following result.

Corollary 3.11**.**

Assume that $\mathbf{M}\in\mathcal{S}(k,\beta,L)$. Under Assumptions 3.1 and 3.6, the choice $N=\lfloor(dTC(\beta,L)/(\|\Sigma_{\varepsilon}\|_{{\rm op}}k))^{1/(2\beta+1)}\rfloor$ ensures

[TABLE]

with probability larger than $1-2e^{-s}$, where $\mathcal{C}>0$ is some constant depending on $\lambda$, $\mathfrak{c}$ and $\mathfrak{K}_{\varepsilon}$.

However, in practice, $\beta$ is not known - nor the rank $k$. This problem is tackled in the next section.

4. Model selection

Assume that we have many possible matrices $\mathbf{\Lambda}_{\tau}$, for $\tau\in\mathcal{T}\subset\{1,\dots,T\}$ and for each $\tau$, many possible $\mathcal{S}_{\tau,k}$ for different possible ranks $k\in\mathcal{K}\subset\{1,\dots,d\wedge T\}$.

Consider $s\in\mathbb{R}_{+}$ and the penalized estimator $\widehat{\mathbf{M}}_{s}=\widehat{\mathbf{M}}_{\mathcal{S}_{\widehat{\tau}_{s},\widehat{k}_{s}}}$ with

[TABLE]

where

[TABLE]

Theorem 4.1**.**

Under Assumptions 3.1 and 3.6, for every $\lambda\in]0,1[$,

[TABLE]

with probability larger than $1-2e^{-s}$.

Remark 4.2**.**

The reader might feel uncomfortable with the fact that the model selection procedure leads to $\hat{k}_{s}$ and $\hat{\tau}_{s}$ that depend on the prescribed confidence level $s$. Note that if $k=k_{0}$ is known, that is $\mathcal{K}=\{k_{0}\}$, then it is clear from the definition that $\hat{\tau}_{s}$ actually does not depend on $s$.

As an application, assume that $\mathbf{M}\in\mathcal{S}(k,\beta,L)$ where $k$ is known, but $\beta$ is unknown. Then the model selection procedure is feasible as it does not depend on $\beta$, and it satisfies exactly the same rate as $\hat{\mathbf{M}}_{\mathcal{S}(k,\beta,L)}$ in Corollary 3.11.

5. Proofs

5.1. Additional notations

Let us first introduce a few additional notations.

First, for the sake of shortness, we introduce the estimation risk $R$ and the empirical risk $r$. These notations also make clear the fact that our estimator can be seen as an empirical risk minimizer.

[TABLE]

Let \Delta(\mathcal{S})=\{\mathbf{A}-\mathbf{B}\textrm{ ; }\mathbf{A},\mathbf{B}\in\mathcal{S}\}.

For any $\mathbf{A}\in\mathcal{M}_{d,T}(\mathbb{C})$, the spectral radius of $\mathbf{A}$ is given by

[TABLE]

Note that $\|\mathbf{A}\|_{\textrm{op}}^{2}=\rho(\mathbf{A}\mathbf{A}^{*})=\rho(\mathbf{A}^{*}\mathbf{A})$.

For any subset $\mathcal{K}$ of $\mathcal{M}_{d,T}(\mathbb{C})$,

[TABLE]

and

[TABLE]

5.2. Some lemmas

Let us now state the key lemmas for the proof of our results. The first one will be used to estimate how far from the minimizer of $R$ is the minimizer of $r$.

Lemma 5.1**.**

For any $\mathbf{A}\in\mathcal{M}_{d,T}(\mathbb{R})$,

[TABLE]

Moreover, for every $\lambda\in]0,1[$,

[TABLE]

and

[TABLE]

Proof.

Consider $\mathbf{A}\in\mathcal{M}_{d,T}(\mathbb{R})$. First of all,

[TABLE]

Then, for any $\lambda\in]0,1[$,

[TABLE]

On the one hand, by Equation (10) together with the classic inequality $2ab\leqslant a^{2}+b^{2}$ for every $a,b\in\mathbb{R}$,

[TABLE]

So, Inequality (8) it true.

On the other hand, by Equation (10) together with the classic inequality $-2ab\leqslant a^{2}+b^{2}$ for every $a,b\in\mathbb{R}$,

[TABLE]

So, Inequality (9) it true. ∎

In the proof of the theorems, $\mathbf{A}$ will be replaced by an estimator of $\mathbf{M}$ that will be data dependent. Thus, it is now crucial to obtain uniform bounds on the scalar product in Lemma 5.1. In machine learning theory, concentration inequalities are the standard tools to derive such a uniform bound, see [6] for a comprehensive introduction to concentration inequalities for independent observations, and their applications to statistics. Some inequalities for time series can be found for example in [19], and were applied to machine learning in [3]. Here, we require more specifically a concentration inequality on random matrices. Such inequalities can be found in  [49, 52]. We will actually use the following result (Theorem 5.39 and Remark 5.40.2 from [52]). As the proof can be found in [52], we don’t reproduce it here.

Proposition 5.2**.**

Under Assumption 3.1, there exists a deterministic constant $\mathfrak{m}>1$, not depending on $\varepsilon$, $d$ and $T$, such that for every $s\in\mathbb{R}_{+}$,

[TABLE]

with probability larger than $1-2e^{-s}$, where $\widetilde{\mathfrak{K}}_{\varepsilon}:=\mathfrak{K}_{\varepsilon}^{2}\vee\mathfrak{K}_{\varepsilon}^{4}$.

We are now in position to provide a uniform bound on the scalar product in Lemma 5.1.

Lemma 5.3**.**

Under Assumption 3.1, there exists a constant $\mathfrak{c}>1$, not depending on $\varepsilon$, $d$ and $T$, such that for every $s\in\mathbb{R}_{+}$ and $\mathcal{K}\subset\mathcal{M}_{d,T}(\mathbb{R})$,

[TABLE]

with probability larger than $1-2e^{-s}$.

Proof.

Consider a subset $\mathcal{K}$ of $\mathcal{M}_{d,T}(\mathbb{R})$ and $s\in\mathbb{R}_{+}$. Let $\sigma_{1}(\varepsilon)\geqslant\dots\geqslant\sigma_{d}(\varepsilon)$ be the singular values of $\varepsilon$. On the one hand, consider a matrix $\mathbf{A}\in\mathcal{K}^{1}$ with singular values $\sigma_{1}(\mathbf{A})\geqslant\dots\geqslant\sigma_{\textrm{rk}(\mathcal{K}^{1})}(\mathbf{A})$. By Cauchy-Schwarz’s inequality:

[TABLE]

Then,

[TABLE]

On the other hand, consider

[TABLE]

Then,

[TABLE]

In particular,

[TABLE]

with $\mathfrak{c}=5\mathfrak{m}/2$. Therefore, by Inequalities (11) and (12) together with Proposition 5.2,

[TABLE]

with probability larger than $1-2e^{-s}$. ∎

5.3. Proof of Theorem 3.2

Consider $\lambda\in]0,1[$ and $s\in\mathbb{R}_{+}$. By applying the Inequalities (8) and (9) of Lemma 5.1 successively:

[TABLE]

By Lemma 5.3:

[TABLE]

with probability larger than $1-2e^{-s}$.

5.4. Proof of Lemma 3.4

First of all,

[TABLE]

Then, since the matrix $\Sigma_{\varepsilon\mathbf{C}}$ is Hermitian,

[TABLE]

5.5. Proof of Lemma 3.5

Since $\widetilde{\varepsilon}=\varepsilon\mathbf{\Lambda}^{+}$, $\Sigma_{\varepsilon}=\sigma^{2}\mathbf{I}_{T}$ and $\mathbf{\Lambda}\mathbf{\Lambda}^{*}=c(\tau,T)\mathbf{I}_{\tau}$,

[TABLE]

Then, for any $x$ with $\|x\|=1$,

[TABLE]

Moreover,

[TABLE]

Therefore,

[TABLE]

and finally, $\widetilde{\mathfrak{K}}_{\widetilde{\varepsilon}}=\widetilde{\mathfrak{K}}_{\varepsilon}$.

5.6. Proof of Theorem 4.1

For short, let us denote

[TABLE]

Consider $\lambda\in]0,1[$. On the one hand, by applying the Inequalities (8) and (9) of Lemma 5.1 successively:

[TABLE]

where

[TABLE]

On the other hand, consider $(\tau,k)\in\mathcal{T}\times\mathcal{K}$. Since $\widetilde{\varepsilon}=\varepsilon\mathbf{\Lambda}_{\tau}^{+}$,

[TABLE]

As in the proof of Proposition 5.3, by Lemma 5.2 and since

[TABLE]

with probability larger than $1-2e^{-u}$,

[TABLE]

Take $u=s+\tau+k$, we obtain that with probability at least $1-2e^{-s-\tau-k}$,

[TABLE]

Then, by a union bound,

[TABLE]

Together with Inequality (5.6), this gives, with probability at least $1-2e^{-s}$,

[TABLE]

Finally, follow the proof of Corollary 3.7 to obtain, on the same event with probability at least $1-2e^{-s}$, for any $\tau$ and $k$,

[TABLE]

Plugging this into (14) gives, with probability at least $1-2e^{-s}$,

[TABLE]

Finally, note that $k\leqslant d$ so

[TABLE]

and so, with probability at least $1-2e^{-s}$,

[TABLE]

This ends the proof.