Blind source separation of tensor-valued time series

TL;DR

This paper extends classic blind source separation methods to tensor-valued time series, providing more effective tools for analyzing high-dimensional data like fMRI and video.

Contribution

It introduces tensorial generalizations of SOBI, FOBI, and JADE, demonstrating their theoretical properties and superior performance over existing methods.

Findings

Tensorial methods outperform vector-based approaches in various settings.

The methods are Fisher consistent and orthogonal equivariant.

Applications to fMRI and video data successfully extract relevant information.

Abstract

The blind source separation model for multivariate time series generally assumes that the observed series is a linear transformation of an unobserved series with temporally uncorrelated or independent components. Given the observations, the objective is to find a linear transformation that recovers the latent series. Several methods for accomplishing this exist and three particular ones are the classic SOBI and the recently proposed generalized FOBI (gFOBI) and generalized JADE (gJADE), each based on the use of joint lagged moments. In this paper we generalize the methodologies behind these algorithms for tensor-valued time series. We assume that our data consists of a tensor observed at each time point and that the observations are linear transformations of latent tensors we wish to estimate. The tensorial generalizations are shown to have particularly elegant forms and we show that each of them is Fisher consistent and orthogonal equivariant. Comparing the new methods with the original ones in various settings shows that the tensorial extensions are superior to both their vector-valued counterparts and to two existing tensorial dimension reduction methods for i.i.d. data. Finally, applications to fMRI-data and video processing show that the methods are capable of extracting relevant information from noisy high-dimensional data.

Equations

$${\mathbf{x}}_t=\mathbf{\Omega }{\mathbf{z}}_t,t=0,\pm 1,\pm 2,\dots ,$$

$${\left(\mathbb{X}{\odot }_m\mathbf{A}\right)}_{i_1\cdots i_r}=\sum _{j_m}{x}_{i_1\cdots i_{m-1}{j}_m{i}_{m+1}\cdots i_r}{a}_{i_m{j}_m},$$

$${\left(\mathbb{X}{\odot }_{-m}\mathbb{Y}\right)}_{kl}=\sum _{i_1,\dots ,i_{m-1},i_{m+1},\dots ,i_r}{x}_{i_1\cdots i_{m-1}k{i}_{m+1}\cdots i_r}{y}_{i_1\cdots i_{m-1}l{i}_{m+1}\cdots i_r}.$$

$${\mathbf{x}}_t^{st}=\mathbf{U}{\mathbf{z}}_t,$$

$${\mathbf{B}}_{\tau ij}({\mathbf{x}}_t)$$

$${\mathbf{\Sigma }}_{\tau }({\mathbf{x}}_t)$$

$${\mathbf{B}}_{\tau }({\mathbf{x}}_t)$$

$${\mathbf{C}}_{\tau ij}({\mathbf{x}}_t)$$

$$\underset{\mathbf{U}{\mathbf{U}}^T=\mathbf{I}}{\text{mbox}max}\sum _{\mathbf{M}\in \mathcal{M}}\Vert \text{mbox}diag({\mathbf{U}}^T\mathbf{M}\mathbf{U}){\Vert }^2,$$

$$\mathbf{\Gamma }({\mathbf{x}}_t):={\mathbf{U}}^T({\mathbf{x}}_t)({\mathbf{\Sigma }}_0({\mathbf{x}}_t){)}^{-1/2},$$

$${\mathbb{X}}_t={\mathbb{Z}}_t{\odot }_1{\mathbf{\Omega }}_1\cdots {\odot }_r{\mathbf{\Omega }}_r,t=0,\pm 1,\pm 2,\dots ,$$

$${\text{x}}_t=({\mathbf{\Omega }}_r^{\ast }⊚\cdots ⊚{\mathbf{\Omega }}_1^{\ast }){\text{z}}_t,$$

$$vec({\mathbb{X}}_t)=({\mathbf{\Omega }}_r\otimes \cdots \otimes {\mathbf{\Omega }}_1)vec({\mathbb{Z}}_t),$$

$${\mathbf{\Sigma }}_0^m({\mathbb{X}}_t):=\frac{1}{{\rho }_m}\text{mbox}E\left[{\mathbf{X}}_t^{(m)}({\mathbf{X}}_t^{(m)}{)}^T\right],$$

$${\mathbb{X}}_t^{st}:={\mathbb{X}}_t{\odot }_1{\left({\mathbf{\Sigma }}_0^1({\mathbb{X}}_t)\right)}^{-1/2}\cdots {\odot }_r{\left({\mathbf{\Sigma }}_0^r({\mathbb{X}}_t)\right)}^{-1/2},$$

$${\mathbb{X}}_t^{st}\propto {\mathbb{Z}}_t{\odot }_1{\mathbf{U}}_1\cdots {\odot }_r{\mathbf{U}}_r,$$

$${\mathbf{B}}_{{\tau }_1{\tau }_2{\tau }_3{\tau }_{4}ij}^m({\mathbb{X}}_t)$$

$${\mathbf{\Sigma }}_{\tau }^m({\mathbb{X}}_t):=$$

$${\mathbf{B}}_{\tau }^m({\mathbb{X}}_t):=$$

$${\mathbf{C}}_{\tau ij}^m({\mathbb{X}}_t):=$$

$$-{\mathbf{\Sigma }}_0^m({\mathbb{X}}_t)({\mathbf{E}}^{ij}+{\mathbf{E}}^{ji}+\mathbf{I}){\mathbf{\Sigma }}_0^m({\mathbb{X}}_t{)}^T,$$

$${\mathbf{\Gamma }}^m({\mathbb{X}}_t)={\mathbf{U}}_m^T({\mathbb{X}}_t)({\mathbf{\Sigma }}_0^m({\mathbb{X}}_t){)}^{-1/2},$$

$$MDI(\hat{\mathbf{\Gamma }})=\frac{1}{\sqrt{p-1}}\underset{\mathbf{C}\in \mathcal{C}}{\inf }\Vert \mathbf{C}\hat{\mathbf{\Gamma }}({\mathbf{A}}_3\otimes {\mathbf{A}}_2\otimes {\mathbf{A}}_1)-\mathbf{I}\Vert ,$$

$${\mathbf{B}}_{{\tau }_1{\tau }_2{\tau }_3{\tau }_{4}ij}^m({\mathbb{X}}_t^{st})=d^4\sum _{k,l}{u}_{ik}^{(m)}{u}_{jl}^{(m)}\cdot {\text{U}}_m{\mathbf{B}}_{{\tau }_1{\tau }_2{\tau }_3{\tau }_{4}kl}^m({\mathbb{Z}}_t){\mathbf{U}}_m^T,$$

$${\mathbf{C}}_{\tau ij}^m({\mathbb{X}}_t^{st})$$

$$+{\mathbf{\Sigma }}_0^m({\mathbb{Z}}_t)[\sum _{k,l}{u}_{ik}^{(m)}{u}_{jl}^{(m)}\cdot ({\mathbf{E}}^{kl}+{\mathbf{E}}^{lk}+\mathbf{I})$$

$$-{\mathbf{U}}_m^T({\mathbf{E}}^{ij}+{\mathbf{E}}^{ji}+\mathbf{I}){\mathbf{U}}_m]{\mathbf{\Sigma }}_0^m({\mathbb{Z}}_t)){\mathbf{U}}_m^T,$$

Full text

††journal: Signal Processing

1 Introduction

1.1 Blind source separation and time series

In the classical blind source separation (BSS) model one assumes that the observed random vectors $\boldsymbol{x}_{i}$, $i=1,\ldots,n$, are linear transformations of some latent vectors of interest, $\boldsymbol{x}_{i}=\boldsymbol{\Omega}\boldsymbol{z}_{i}$, where $\boldsymbol{\Omega}\in\mathbb{R}^{p\times p}$ is a full rank mixing matrix. Coupling the model with different sets of assumptions on $\boldsymbol{z}_{i}$ gives various well-known models: (i) assuming that $\boldsymbol{z}_{i}$ are i.i.d. and have mutually independent components yields the independent component (IC) model [see e.g. Comon and Jutten, 2010]; (ii) assuming $\boldsymbol{z}_{i}$ are i.i.d. and have spherical distribution yields an elliptical model for $\boldsymbol{x}_{i}$ [Oja, 2010] and (iii) as a special case of both previous, assuming that $\boldsymbol{z}_{i}$ are i.i.d. and have standard Gaussian distribution yields the general multivariate Gaussian distribution for $\boldsymbol{x}_{i}$.

In the context of time series it is natural to incorporate the time dependency of the data into the model structure and a commonly used BSS model [Belouchrani et al., 1997] assumes that the observed series $\boldsymbol{x}_{t}$ is generated as

[TABLE]

where the latent time series $\boldsymbol{z}_{t}$ satisfies the following three assumptions.

(V1)****.

$\mbox{E}[\boldsymbol{z}_{t}]=\boldsymbol{0}$**

(V2)****.

$\mbox{Cov}[\boldsymbol{z}_{t}]=\boldsymbol{I}$**

(V3)****.

$\mbox{E}[\boldsymbol{z}_{t}\boldsymbol{z}_{t+\tau}^{T}]=\mbox{E}[\boldsymbol{z}_{t+\tau}\boldsymbol{z}_{t}^{T}]=\boldsymbol{D}_{\tau}$* is diagonal for all $\tau=1,2,\ldots$.*

Without loss of generality, (V1) implies that the observed series is centered and (V2) fixes the scales of the columns of $\boldsymbol{\Omega}$. After these two assumptions we can still freely change the signs and order of the elements of $\boldsymbol{z}_{t}$ and the corresponding columns of $\boldsymbol{\Omega}$ without altering the overall model. Thus the order and the signs of the latent series are unidentifiable which, however, is usually not a problem in practice. (V1)–(V3) together also imply that the time series $\boldsymbol{z}_{t}$ and $\boldsymbol{x}_{t}$ are weak second-order stationary and $\boldsymbol{x}_{t}$ satisfies $\mbox{E}[\boldsymbol{x}_{t}]=\boldsymbol{0}$ and $\mbox{Cov}[\boldsymbol{x}_{t}]=\boldsymbol{\Omega}\boldsymbol{\Omega}^{T}$.

Time series BSS models such as this have gained in popularity in recent years as general multivariate time series models are demanding in both theory and computation. On the other hand, assuming a BSS model allows the use of well-established univariate time series methods for each of the estimated latent components separately. For some recent contributions see, for example, Oja et al. [2000], Wu and Yu [2005], Chen et al. [2007], Broda and Paolella [2009], Lu et al. [2009], Chen et al. [2010], García-Ferrer et al. [2012].

1.2 Tensor-valued methods for time series

In the models discussed above at each time point a $p$-variate vector is observed. Modern data structures are however often more complex. For example, in many applications at each time point data might be observed which is better represented by a tensor. Such applications are for instance spatio-temporal data where at each time point usually a matrix is obtained or fMRI (functional Magnetic Resonance Imaging) data where for each time point a 3-dimensional tensor is recorded. But also video clip data can be seen as a time series where each frame is a matrix- or tensor-valued observation, depending on the number of colors used.

The most common approach to analyzing such data is to convert, following some convention, the tensor into a large vector and then apply standard multivariate methods for vector-valued data. Besides the often practical problem that the vectorized data might be of quite high dimension also information gets lost in this process. As for example Werner et al. [2008] point out, after vectorizing the tensors the resulting vectors have a Kronecker structure. Ignoring this structure then means that a much larger number of parameters needs to be estimated.

For i.i.d. data this has recently led to extensive research where methods either model the Kronecker structure or work directly with the tensors. For some recent work on structured multivariate estimation see for example Werner et al. [2008], Srivastava et al. [2008], Wiesel [2012], Greenewald and Hero [2015] and references therein. For some recent contributions for i.i.d. tensor methods see for example Li et al. [2010], Zeng and Zhong [2013], Ding and Cook [2015], Zhong et al. [2015] and references therein.

Also independent component analysis (ICA) has already been considered in the context of tensors, some early works being Beckmann and Smith [2005], Zhang et al. [2008], Vasilescu and Terzopoulos [2005]. A fully tensorial model-based approach was however only recently developed in Virta et al. [2016a, b] where tensorial versions of the well-known ICA methods Fourth order blind identification (FOBI) [Cardoso, 1989] and Joint approximate diagonalization of eigen-matrices (JADE) [Cardoso and Souloumiac, 1993] were introduced.

Methods for tensor-valued time series seem however not to have had much attention yet although the first steps are for example Walden and Serroukh [2002], Rogers et al. [2013], Bahadori et al. [2014]. But to the best of our knowledge no tensorial BSS-methods for dependent data have been considered so far. To fill this gap, in this paper we propose tensor extensions for three BSS methods meant for multivariate time series. The first method is called the Second order blind identification (SOBI) [Belouchrani et al., 1997] and is based on using second-order information in the form of autocovariance matrices to separate the hidden source series. As such it is best suited for multivariate linear processes and may not work with models having trivial autocovariances, for example with stochastic volatility models such as GARCH. The recently proposed second and third methods, generalized FOBI (gFOBI) and generalized JADE (gJADE), are the exact opposite and operate on the joint fourth-order moments of the component series, see Matilainen et al. [2015]. Thus the successful use of either requires non-trivial higher moments, ruling out for example the standard ARMA models. The tensorial extensions of these three methods are respectively called TSOBI, TgFOBI and TgJADE and are discussed in Section 4.

1.3 Structure of the paper

In Section 2 we introduce the used notation and define various concepts of multilinear algebra we need to operate the tensor observations. Although fairly easy to grasp and use, after we define the $m$-flattening of a tensor most tensor operations can be carried out conveniently in a matrix form. Section 3 reviews the theory of SOBI, gFOBI and gJADE and prepares the ground for their tensor versions in Section 4 where the corresponding theory and algorithms are discussed. In Section 5 we first use simulations to compare the presented methods with their vector-valued counterparts for vectorized data and then use the methods to process simulated fMRI-data and a color video. In both applications the proposed methodology is shown to extract the key elements of the signals in compressed form. In Section 6 we finally conclude with some discussion and the proofs are gathered in A.

2 Notation and tensor algebra

2.1 Notation in general

Throughout the paper scalars are denoted by lower-case letters, $a,b,c,\ldots$, vectors by lower-case boldface letters, $\boldsymbol{a},\boldsymbol{b},\boldsymbol{c},\ldots$, matrices by capital boldface letters, $\boldsymbol{A},\boldsymbol{B},\boldsymbol{C},\ldots$, and tensors of general order by capital blackboard letters, $\mathbb{A},\mathbb{B},\mathbb{C},\ldots$ (note that $\mathbb{R}$ still means the real line). The same convention on fonts is followed with random elements, but instead using the letters from the end of the alphabet, $x,y,z,\boldsymbol{x},\boldsymbol{y},\boldsymbol{z}$, etc.

2.2 Regarding matrices

The standard basis vectors of $\mathbb{R}^{p}$ are denoted by $\boldsymbol{e}_{i}$, $i=1,\ldots,p$, and using them we can construct the matrices $\boldsymbol{E}^{ij}:=\boldsymbol{e}_{i}\boldsymbol{e}_{j}^{T}$, the only non-zero element of $\boldsymbol{E}^{ij}$ being the single one as the element $(i,j)$. We further make use of the following sets of $p\times p$ matrices: $\mathcal{P}$, the set of all matrices with a single one in each row and column and rest of the entries zero; $\mathcal{J}$, the set of all diagonal matrices with the diagonal entries equal to $\pm 1$; $\mathcal{D}$, the set of all diagonal matrices with positive diagonal elements and $\mathcal{C}$, the set of all matrices $\boldsymbol{C}=\boldsymbol{PJD}$ where $\boldsymbol{P}\in\mathcal{P}$, $\boldsymbol{J}\in\mathcal{J}$ and $\boldsymbol{D}\in\mathcal{D}$. The sets $\mathcal{P},\mathcal{J}$ and $\mathcal{D}$ then respectively correspond to the sets of permutation matrices, heterogeneous sign-change matrices and heterogeneous scaling matrices. Finally, $\|\cdot\|$ is the Frobenius norm and by the equivalence $\boldsymbol{A}\equiv\boldsymbol{B}$ we mean that $\boldsymbol{A}=\boldsymbol{PJB}$ for some $\boldsymbol{P}\in\mathcal{P}$ and $\boldsymbol{J}\in\mathcal{J}$.

2.3 Regarding tensors

To manipulate tensors we next provide some basic tools of multilinear algebra, see also De Lathauwer et al. [2000]. By a tensor-valued time series $\mathbb{X}_{t}$ we mean the set $\{\mathbb{X}_{t}\}_{t=-T}^{T}$ of realisations of a tensor-valued stochastic process $\mathbb{X}_{t}\in\mathbb{R}^{p_{1}\times\cdots\times p_{r}}$ on some fixed set of time points $t=-T,\ldots,T$. That is, for each time point we observe a tensor of the same size, something akin to the frames of a video, see Figure 1 for a visual representation.

A useful way of understanding a tensor as a collection of blocks of lower order is provided by the $m$-mode vectors (or fibers), which extend the notion of the columns and rows of a matrix into higher dimensions. The $m$-mode vectors of a tensor $\mathbb{X}\in\mathbb{R}^{p_{1}\times\cdots\times p_{r}}$ are obtained by fixing the values of all other indices but the $m$th and letting the $m$th index vary over its range. The columns and rows of a matrix are thus also its $1$-mode and $2$-mode vectors, respectively. Each $m$-mode vector is of length $p_{m}$ and the total number of $m$-mode vectors in $\mathbb{X}$ is $\rho_{m}:=\prod_{i\neq m}^{r}p_{i}$. Stacking all of them horizontally in a pre-specified order into a matrix yields what we call the $m$-flattening of a tensor, denoted by $\boldsymbol{X}^{(m)}\in\mathbb{R}^{p_{m}\times\rho_{m}}$. As is later seen the ordering is not important as long as it is consistent and one option is e.g. the cyclical ordering discussed in De Lathauwer et al. [2000].

We next define two multiplication operations, the first one being the $m$th mode linear transformation $(\mathbb{X}\odot_{m}\boldsymbol{A})\in\mathbb{R}^{p_{1}\times\cdots\times p_{r}}$ defined element-wise as

[TABLE]

where $\mathbb{X}\in\mathbb{R}^{p_{1}\times\cdots\times p_{r}}$ and $\boldsymbol{A}\in\mathbb{R}^{p_{m}\times p_{m}}$. The second operation, $(\mathbb{X}\odot_{-m}\mathbb{Y})\in\mathbb{R}^{p_{m}\times p_{m}}$, takes two tensors of the same size, $\mathbb{X},\mathbb{Y}\in\mathbb{R}^{p_{1}\times\cdots\times p_{r}}$ and acts element-wise as

[TABLE]

Nice interpretations for the previous operations are given by noticing that $(\mathbb{X}\odot_{m}\boldsymbol{A})^{(m)}=\boldsymbol{A}\boldsymbol{X}^{(m)}$ and that $\mathbb{X}\odot_{-m}\mathbb{X}=\boldsymbol{X}^{(m)}(\boldsymbol{X}^{(m)})^{T}$. That is, $\mathbb{X}\odot_{m}\boldsymbol{A}$ applies the linear transformation given by $\boldsymbol{A}$ individually to all $m$-mode vectors of $\mathbb{X}$ and $\mathbb{X}\odot_{-m}\mathbb{X}$ is simply the sum of the outer products of all $m$-mode vectors of $\mathbb{X}$ with themselves, the order we stacked them into $\boldsymbol{X}^{(m)}$ clearly not mattering.

Although the above tensor notation is quite straightforward, the familiarity of matrix notation still outweighs its usefulness and we thus choose to work with the $m$-flattened tensors whenever possible. The ability to express everything with matrices also makes the implementation of the algorithms easier.

3 Blind source separation of multivariate time series

3.1 The model and the standardized time series

In this section we review the theory behind three moment-based blind source separation methods for time series, SOBI, gFOBI and gJADE, which are then in Section 4 extended for tensor-valued time series.

All three methods assume the model (1) along with (V1) and (V2) but each of them has additional assumptions that guarantee the successful identification of the latent time series with the particular method. These extra assumptions will be elaborated more on the next section discussing the estimation of the rotation.

Interestingly, the model and (V1) and (V2) alone already yield the first step in the estimation of $\boldsymbol{z}_{t}$. This step is given by the whitening of the data, $\boldsymbol{x}_{t}\mapsto\boldsymbol{x}_{t}^{st}:=(\boldsymbol{\Sigma}_{0}(\boldsymbol{x}_{t}))^{-1/2}\boldsymbol{x}_{t}$, where $\boldsymbol{\Sigma}_{0}(\boldsymbol{x}_{t})$ is the covariance matrix of the stationary series $\boldsymbol{x}_{t}$. It is easy to see (e.g. by considering the singular value decomposition of $\boldsymbol{\Omega}$) that the whitened series satisfies

[TABLE]

for some orthogonal matrix $\boldsymbol{U}\in\mathbb{R}^{p\times p}$. The problem is then reduced from that of estimating the $p^{2}$ parameters of $\boldsymbol{\Omega}^{-1}$ into that of estimating the $(1/2)p(p-1)$ parameters of $\boldsymbol{U}$, the standardization thus solving half of the problem in terms of number of parameters to estimate.

3.2 Estimating the rotation

The unknown rotation $\boldsymbol{U}$ is estimated using one of three families of matrix-valued functionals. But prior to defining them we first define another set of functionals that is used in constructing one of the families. $\mathcal{T}$ being a pre-specified set of lags, we define

[TABLE]

where $\tau\in\mathcal{T}$ and $i,j=1,\ldots,p$. The matrices $\boldsymbol{B}_{\tau ij}(\boldsymbol{x}^{st}_{t})$ then contain specific sums of joint lagged fourth-order moments of the series $\boldsymbol{x}^{st}_{t}$ and, concerning the current estimation problem, unless the underlying series has trivial fourth-order moments they contain information on the unknown rotation.

The three families of interest are then defined as follows.

[TABLE]

where $\tau\in\mathcal{T}$, $i,j=1,\ldots,p$ and $\delta_{ij}$ is the Kronecker delta. Consider then the following three additional assumptions.

(V4)****.

The matrix $\boldsymbol{\Sigma}_{\tau}(\boldsymbol{z}_{t})$ is diagonal for all $\tau\in\mathcal{T}$.

(V5)****.

The matrix $\boldsymbol{B}_{\tau}(\boldsymbol{z}_{t})$ is diagonal for all $\tau\in\mathcal{T}$.

(V6)****.

The matrix $\boldsymbol{C}_{\tau ij}(\boldsymbol{z}_{t})$ is diagonal for all $\tau\in\mathcal{T}$, $i,j=1,\ldots,p$.

If we then assume the model (1), (V1) and (V2) and additionally either (V4), (V5) or (V6) it can be shown that the matrix $\boldsymbol{U}^{T}$ diagonalizes respectively the matrices $\boldsymbol{\Sigma}_{\tau}(\boldsymbol{x}^{st}_{t})$, $\boldsymbol{B}_{\tau}(\boldsymbol{x}^{st}_{t})$ or $\boldsymbol{C}_{\tau ij}(\boldsymbol{x}^{st}_{t})$ for all $\tau\in\mathcal{T}$, $i,j=1,\ldots,p$. Technically the eigendecomposition of any one of them could then under the suitable assumptions be used to estimate the needed rotation. However, for the rotation to be unique the eigenvalues of the considered matrix need to be distinct and this is not necessarily the case. One solution is, instead of diagonalizing just one of the matrices, to approximately diagonalize several of them using the so-called joint approximate diagonalization. A often used technique for performing the joint diagonalization of a chosen set of matrices $\mathcal{M}$ is given by the optimization problem

[TABLE]

that is, we find the orthogonal matrix that makes all matrices in the set as diagonal as possible, as measured by the sum of squares of the resulting diagonal values. It can be shown that this is equivalent to minimizing the sum of squared off-diagonal values. Choosing $\mathcal{M}=\{\boldsymbol{\Sigma}_{\tau}(\boldsymbol{x}^{st}_{t})\}_{\tau\in\mathcal{T}}$, $\mathcal{M}=\{\boldsymbol{B}_{\tau}(\boldsymbol{x}^{st}_{t})\}_{\tau\in\mathcal{T}}$ or $\mathcal{M}=\{\boldsymbol{C}_{\tau ij}(\boldsymbol{x}^{st}_{t})\}_{\tau\in\mathcal{T},i,j=1,\ldots,p}$ gives then SOBI, gFOBI and gJADE, respectively. For further details on the joint diagonalization, see Belouchrani et al. [1997], Illner et al. [2015].

Intuitively, diagonalizing a set of matrices means that no single matrix needs to carry the information to separate all $p$ components as long as the information exists in at least one of the matrices in the set. This reasoning is encoded in the following three assumptions that guarantee that the correct solution is actually found in SOBI, gFOBI and gJADE, respectively.

(V7)****.

For all pairs $i\neq j$ there exists $\tau\in\mathcal{T}$ such that the $i$th and $j$th diagonal elements of $\boldsymbol{\Sigma}_{\tau}(\boldsymbol{z}_{t})$ are distinct.

(V8)****.

For all pairs $i\neq j$ there exists $\tau\in\mathcal{T}$ such that the $i$th and $j$th diagonal elements of $\boldsymbol{B}_{\tau}(\boldsymbol{z}_{t})$ are distinct.

(V9)****.

For $p-1$ components $i$ there exists $\tau\in\mathcal{T}$ such that $\boldsymbol{C}_{\tau ii}(\boldsymbol{z}_{t})\neq\boldsymbol{0}$.

The nature of the above assumptions shows that SOBI is targeted to a different set of distributions for $\boldsymbol{z}$ than gFOBI and gJADE; SOBI uses second order information in the form of autocovariances to separate the latent time series while gFOBI and gJADE implicitly assume no second order information exists and move directly to the use fourth-order moments.

To put everything together, all algorithms require the model (1) along with (V1) and (V2) and additionally SOBI assumes (V4) and (V7); gFOBI assumes (V5) and (V8) and gJADE assumes (V6) and (V9). Note that the standard assumption (V3) implies the weaker one (V4) and can replace it if needed. Also, since (V5) and (V6) have little meaning in applications sometimes it might be more practical to assume the following stronger condition that implies all of (V3), (V4), (V5) and (V6).

(V10)****.

The component series of $\boldsymbol{z}_{t}$ are independent.

Having estimated the rotation with one of the methods we define the related unmixing matrix functional as

[TABLE]

where $\boldsymbol{U}^{T}(\boldsymbol{x}_{t})$ is the diagonalizing rotation estimated from the series $\boldsymbol{x}_{t}$ notated as a functional of the data. The estimated components are then obtained via the transformation $\boldsymbol{x}_{t}\mapsto\boldsymbol{\Gamma}(\boldsymbol{x}_{t})\boldsymbol{x}_{t}$. One important property of SOBI, gFOBI and gJADE is that they are all affine equivariant, that is, the related unmixing matrix functionals satisfy $\boldsymbol{\Gamma}(\boldsymbol{A}\boldsymbol{x}_{t})\equiv\boldsymbol{\Gamma}(\boldsymbol{x}_{t})\boldsymbol{A}^{-1}$ for all non-singular $\boldsymbol{A}$, even for $\boldsymbol{x}_{t}$ outside the model (1). Thus the choice of the coordinate system does not affect the resulting estimated component series and the three functionals can be seen as examples of invariant coordinate system (ICS) functionals, see Tyler et al. [2009].

4 Tensor blind source separation

4.1 The model and the standardized time series

Using the linear transformation operator $\odot_{m}$ we define the tensor extension of the model (1) in a manner similar to the ICA model extensions in Virta et al. [2016a, b]. This yields the following tensor blind source separation model for the tensor-valued time series $\mathbb{X}_{t}$.

[TABLE]

where $\boldsymbol{\Omega}_{m}\in\mathbb{R}^{p_{m}\times p_{m}}$, $m=1,\ldots,r$, are non-singular $m$-mode mixing matrices and the unobserved tensor-valued time series $\mathbb{Z}_{t}$ satisfies the following three assumptions for all $t=0,\pm 1,\pm 2,\ldots$

(T1)****.

$\mbox{E}[vec(\mathbb{Z}_{t})]=\boldsymbol{0}$.

(T2)****.

$\mbox{Cov}[vec(\mathbb{Z}_{t})]=\boldsymbol{I}$.

(T3)****.

$\mbox{E}[\mathbb{Z}_{t}^{(m)}(\mathbb{Z}_{t+\tau}^{(m)})^{T}]=\mbox{E}[\mathbb{Z}_{t+\tau}^{(m)}(\mathbb{Z}_{t}^{(m)})^{T}]=\boldsymbol{D}^{m}_{\tau}$* are diagonal for all $m=1,\ldots,r$ and for all $\tau=1,2,\ldots$*

Assumption (T1) again without loss of generality requires that the observed time series is centered and Assumption (T2) fixes the relative scales of the columns of the mixing matrices. Note that one can still multiply any of the mixing matrices by a constant and divide one of them by the same constant without changing the model. Assumption (T2) further implies that the covariance matrix of the vectorized series $vec(\mathbb{X}_{t})$ has the Kronecker structure, $\mbox{Cov}(vec(\mathbb{X}_{t}))=(\boldsymbol{\Omega}_{m}\boldsymbol{\Omega}_{m}^{T})\otimes\cdots\otimes(\boldsymbol{\Omega}_{1}\boldsymbol{\Omega}_{1}^{T})$, where $\otimes$ is the Kronecker product. Assumption (T3) is a straight generalization of (V3) and, recalling the definition of $\boldsymbol{Z}^{(m)}_{t}$, in principle says that the $m$-mode vectors of $\mathbb{Z}_{t}$ are “on average stationary”. Like for the vector methods, also now we need additional assumptions to guarantee that each of the methods finds a solution and that the solution is the correct one and these are discussed in the next section.

Curiously, the BSS model (9) can be viewed as a collection of low-multilinear rank models for each observed tensor $\mathbb{Z}_{t}$ with the constraint that the factor matrices are equal for each time point. Still another viewpoint is to treat model (9) as a low-multilinear rank model for the whole tensor-valued time series viewed as a single tensor of order $r+1$, with the time dimension left uncompressed. If we further impose the additional constraint that the mixing matrices be orthogonal the model can be seen as a multilinear singular value decomposition in the same two different ways.

Two interesting analogies between the model (9) and tensor segmentation can also be drawn. In Boussé et al. [2017] it is assumed that in vector-valued BSS either the individual signals or the columns of the mixing matrix (and later both) admit low-rank tensorizations. Consider first the former case and assume that we have a data matrix $\textbf{X}\in\mathbb{R}^{T\times p}$ with the $p$ columns corresponding to the $p$ individual observed signals. Then our working assumption is that the rows of X are actually vectorized tensors while the segmentation conversely assumes that the columns of X are actually vectorized tensors. Consider next the case where the columns of the mixing matrix are actually vectorized tensors and therein a subcase where each of the tensors is of rank one. Then the vector BSS-model (1) can be written as

[TABLE]

where the vectors and matrices are of suitable size, the $j$th columns of the matrices $\boldsymbol{\Omega}^{*}_{m}$, $m=1,\ldots,r$, give the factors of the rank-1 tensorization of the $j$th column of the mixing matrix $\boldsymbol{\Omega}$, $j=1,\ldots,p$ and $\circledcirc$ denotes the Khatri-Rao product (column-wise Kronecker product). Comparing (10) with the vectorization of the model (9),

[TABLE]

shows that the two approaches are linked, both assuming a certain structure on the mixing matrix.

As in the vector case, the first step towards solving (9) is given by a standardization based on second-order moments. Define first the $m$th mode covariance matrices as

[TABLE]

where $m=1,\ldots,r$. The matrix $\boldsymbol{\Sigma}_{0}^{m}(\mathbb{X}_{t})$ is essentially the mean covariance matrix of all $m$-mode vectors of $\mathbb{X}_{t}$ and thus measures the average spatial dependency in the $m$th direction. The first step of the unmixing procedure is then given by a simultaneous standardization from all modes of $\mathbb{X}_{t}$:

[TABLE]

where $\left(\boldsymbol{\Sigma}_{0}^{m}(\mathbb{X}_{t})\right)^{-1/2}$ is the unique symmetric inverse square root of $\boldsymbol{\Sigma}_{0}^{m}(\mathbb{X}_{t})$, $m=1,\ldots,r$. This particular standardization was first considered by Virta et al. [2016a] in the context of independent component analysis and there the authors also gave the following theorem (in the guise of the ICA model).

Theorem 1**.**

Assume that the tensor-valued time series $\mathbb{X}_{t}$ is generated by the model (9) and satisfies Assumptions (T1) and (T2). Then

[TABLE]

for some orthogonal matrices $\boldsymbol{U}_{1}\in\mathbb{R}^{p_{1}\times p_{1}}$, …, $\boldsymbol{U}_{r}\in\mathbb{R}^{p_{r}\times p_{r}}$. The constant of proportionality is $(\prod_{m=1}^{r}p_{m}^{1/2})^{r-1}\|\boldsymbol{\Omega}_{r}\otimes\cdots\otimes\boldsymbol{\Omega}_{1}\|^{1-r}$.

The proof of the above theorem does not depend in any way on the time series nature of $\mathbb{X}_{t}$ and is essentially similar to that of Theorem 5.3.1 in Virta et al. [2016a] and thus omitted here. Comparing the result of Theorem 1 to (2) shows that the former serves as a direct tensor analogy to the latter and again reduces the problem of inverting all matrices $\boldsymbol{\Omega}_{m}$ into the problem of inverting the unknown rotations $\boldsymbol{U}_{m}$.

4.2 Estimating the rotations

For estimating the unknown rotations we use tensor counterparts for the families of functionals (4), (5) and (6), but first we have to define a counterpart for (3). Assuming a fixed mode $m$ and fixed lags $\tau_{1}$, $\tau_{2}$, $\tau_{3}$ and $\tau_{4}$ we define

[TABLE]

where $i,j=1,\ldots,p_{m}$. Using the above we can then define the functionals that allow the estimation of the $r$ rotations. Each of the following is a $p_{m}\times p_{m}$ matrix-valued functional and we again assume that $\mathcal{T}$ is a fixed set of lags.

[TABLE]

where $\tau\in\mathcal{T}$ and $i,j=1,\ldots,p_{m}$. The matrices $\boldsymbol{B}_{\tau}^{m}(\mathbb{X}_{t})$ and $\boldsymbol{C}_{\tau ij}^{m}(\mathbb{X}_{t})$ can be seen as lagged versions of the matrices used in TFOBI and TJADE, see Virta et al. [2016a, b]. Note that the basic form of the matrices $\boldsymbol{C}_{\tau ij}^{m}(\mathbb{X}_{t})$ in (17) differs from the form of their vector counterparts in (6). This modification is needed to obtain the results in Theorem 2 below.

We further introduce the following three assumptions, relating to the matrices (15), (16) and (17), respectively.

(T4)****.

For all $m=1,\ldots,r$ and $\tau\in\mathcal{T}$ the matrix $\boldsymbol{\Sigma}^{m}_{\tau}(\mathbb{Z}_{t})$ is diagonal.

(T5)****.

For all $m=1,\ldots,r$ and $\tau\in\mathcal{T}$ the matrix $\boldsymbol{B}^{m}_{\tau}(\mathbb{Z}_{t})$ is diagonal.

(T6)****.

For all $m=1,\ldots,r$, $i,j=1,\ldots,p_{m}$ and $\tau\in\mathcal{T}$ the matrix $\boldsymbol{C}^{m}_{\tau ij}(\mathbb{Z}_{t})$ is diagonal.

Putting all the previous concepts together, a connection between the functionals (15), (16) and (17) and the estimation of the rotations is stated in the next theorem.

Theorem 2**.**

Assume that the tensor-valued time series $\mathbb{X}_{t}$ is generated by the model (9) and satisfies Assumptions (T1), (T2) and additionally either $i)$ (T4), $ii)$ (T5) or $iii)$ (T6). Then, for a fixed mode $m=1,\ldots,r$, the matrix $\boldsymbol{U}_{m}^{T}$, where $\boldsymbol{U}_{m}$ is the $m$th true rotation in Theorem 1, diagonalizes respectively $i)$ the matrices $\boldsymbol{\Sigma}_{\tau}^{m}(\mathbb{X}^{st}_{t})$, $\tau\in\mathcal{T}$, $ii)$ the matrices $\boldsymbol{B}_{\tau}^{m}(\mathbb{X}^{st}_{t})$, $\tau\in\mathcal{T}$ or $iii)$ the matrices $\boldsymbol{C}_{\tau ij}^{m}(\mathbb{X}^{st}_{t})$, $\tau\in\mathcal{T}$, $i,j=1,\ldots,p_{m}$.

Following the ideas of Section 3, we can then estimate the $m$th rotation $\boldsymbol{U}_{m}^{T}$ by simultaneously diagonalizing, see (7), one of the matrix sets: $\mathcal{M}_{m}=\{\boldsymbol{\Sigma}_{\tau}^{m}(\mathbb{X}^{st}_{t})\}_{\tau\in\mathcal{T}}$, $\mathcal{M}_{m}=\{\boldsymbol{B}_{\tau}^{m}(\mathbb{X}^{st}_{t})\}_{\tau\in\mathcal{T}}$ or $\mathcal{M}_{m}=\{\boldsymbol{C}^{m}_{\tau ij}(\mathbb{X}^{st}_{t})\}_{\tau\in\mathcal{T},i,j=1,\ldots,p}$, yielding TSOBI, TgFOBI and TgJADE, respectively. Note however that Theorem 2 does not guarantee that the estimated diagonalizing rotation is necessarily the correct one (it might be a different basis for the correct eigenspace) and to ensure that we further need the following assumptions which apply to TSOBI, TgFOBI and TgJADE, respectively.

(T7)****.

For all $m=1,\ldots,r$ and for all pairs $i\neq j$ there exists $\tau\in\mathcal{T}$ such that the $i$th and $j$th diagonal elements of $\boldsymbol{\Sigma}^{m}_{\tau}(\mathbb{Z}_{t})$ are distinct.

(T8)****.

For all $m=1,\ldots,r$ and for all pairs $i\neq j$ there exists $\tau\in\mathcal{T}$ such that the $i$th and $j$th diagonal elements of $\boldsymbol{B}^{m}_{\tau}(\mathbb{Z}_{t})$ are distinct.

(T9)****.

For all $m=1,\ldots,r$, for $p_{m}-1$ indices $i$ there exists $\tau\in\mathcal{T}$ such that $\boldsymbol{C}^{m}_{\tau ii}(\mathbb{Z}_{t})\neq\boldsymbol{0}$.

Putting the assumptions again together all three methods assume the model (9) with (T1) and (T2). Additionally, TSOBI requires (T4) and (T7); TgFOBI requires (T5) and (T8) and TgJADE requires (T6) and (T9). Since the assumptions (T4), (T5) and (T6) regarding the working of the methods are again quite impractical it is worth noting that the general assumption (T3) again implies (T4) and all of (T3), (T4), (T5) and (T6) are implied by the following stronger, but more intuitive assumption.

(T10)****.

The component series of $\mathbb{Z}_{t}$ are independent.

This implication can be proven by inspecting the matrices $\boldsymbol{\Sigma}^{m}_{\tau}(\mathbb{Z}_{t})$, $\boldsymbol{B}^{m}_{\tau}(\mathbb{Z}_{t})$ and $\boldsymbol{C}^{m}_{\tau ij}(\mathbb{Z}_{t})$ element-wise under (T10), e.g. in the manner of proof of Theorem 1 in Virta et al. [2016b].

An algorithm for applying the three proposed methods is given next. The choices for the set of matrices $\mathcal{M}_{m}$ leading to the different methods are listed above after Theorem 2.

Remark 1**.**

Like in the original SOBI, the matrices $\boldsymbol{\Sigma}^{m}_{\tau}(\mathbb{X}^{st}_{t})$ are in practice rarely symmetric, thus not satisfying (T4), and in the algorithm one should instead use the symmetrized versions, $(1/2)(\boldsymbol{\Sigma}^{m}_{\tau}(\mathbb{X}^{st}_{t})+\boldsymbol{\Sigma}^{m}_{\tau}(\mathbb{X}^{st}_{t})^{T})$. Whichever we use makes theoretically no difference as both are under the model diagonalized by the same orthogonal matrix. This practice is described also for other SOBI variants e.g. in Miettinen et al. [2014], Ilmonen et al. [2015].

As in the vector case, we can for tensor-valued time series define the $m$-mode unmixing functionals as

[TABLE]

where $\boldsymbol{U}_{m}^{T}(\mathbb{X}_{t})$ is the diagonalizing $m$th rotation estimated from the series $\mathbb{X}_{t}$ notated as a functional of the data. Using the unmixing functionals of all $r$ modes the final solution is then obtained via the transformation $\mathbb{X}_{t}\mapsto\mathbb{X}_{t}\odot_{1}\boldsymbol{\Gamma}^{1}(\mathbb{X}_{t})\cdots\odot_{r}\boldsymbol{\Gamma}^{r}(\mathbb{X}_{t})$. Of the tensor extensions, however, none is affine equivariant in the same sense as their vector counterparts if $r>1$. That is, it is not in general true that $\boldsymbol{\Gamma}^{m}(\mathbb{X}_{t}\odot_{1}\boldsymbol{A}_{1}\cdots\odot_{r}\boldsymbol{A}_{r})\equiv\boldsymbol{\Gamma}^{m}(\mathbb{X}_{t})\boldsymbol{A}_{m}^{-1}$ for all non-singular $\boldsymbol{A}_{1},\ldots,\boldsymbol{A}_{r}$. However, the previous holds when all $\boldsymbol{A}_{1},\ldots,\boldsymbol{A}_{r}$ are orthogonal and all the tensorial extensions are thus orthogonally equivariant. See Virta et al. [2016a] for more discussion and a conjecture that affine equivariance is unobtainable in the general tensor case.

Remark 2**.**

Using no temporal information, that is, choosing $\mathcal{T}=\{0\}$ in either TgFOBI or TgJADE yields respectively the tensorial ICA methods TFOBI and TJADE as introduced in Virta et al. [2016a, b].

Remark 3**.**

In Virta et al. [2016a, b] it is shown that both TFOBI and TJADE have alternative versions for (16) and (17) and such modifications could easily be formulated also for TgFOBI and TgJADE. However, Virta et al. [2016a] also showed using asymptotic properties that for the alternative form of TFOBI to have lower total asymptotic variance than the main form the majority of the components of the latent tensor need to have nearly symmetric distributions. Virta et al. [2016b] on the other hand showed that both versions of TJADE produce almost identical results, the alternative form being computationally much slower. Based on these considerations we thus choose to work with the primary forms in (16) and (17).

5 Examples

5.1 Simulation study

We next compared the performances of the proposed methods to those of their vector-valued counterparts by applying the former to the actual simulated tensor-valued data and the latter to the same data in vectorized form. The proposed methods take advantage of the tensor structure and are thus expected to perform better. The simulations were done with R [R Core Team, 2015] using the packages fGarch [Wuertz et al., 2013], ggplot2 [Wickham, 2009], JADE [Miettinen et al., 2017], stochvol [Kastner, 2016] and tsBSS [Matilainen et al., 2016]. Additionally, our package tensorBSS [Virta et al., 2016c] provides implementations of all the discussed tensor methods.

We considered two settings both consisting of simulated tensor-valued time series $\mathbb{Z}_{t}\in\mathbb{R}^{3\times 2\times 2}$ with the lengths of either $T=1000$, $2000$, $4000$, $8000$, $16000$ or $32000$. In the first setting all component series of $\mathbb{Z}_{t}$ were generated from different ARMA models and in the second setting from various models exhibiting stochastic volatility (SV), see A for detailed information. The tensors $\mathbb{Z}_{t}$ were then mixed as $\mathbb{X}_{t}:=\mathbb{Z}_{t}\odot_{1}\boldsymbol{A}_{1}\odot_{2}\boldsymbol{A}_{2}\odot_{3}\boldsymbol{A}_{3}$ where the square matrices $\boldsymbol{A}_{1}$, $\boldsymbol{A}_{2}$ and $\boldsymbol{A}_{3}$ all either had elements sampled independently from the standard normal distribution or were random orthogonal matrices uniform with respect to the Haar measure.

In addition to the proposed methods also the non-lagged versions of TgFOBI and TgJADE, namely TFOBI and TJADE, along with their vector counterparts FOBI and JADE were included in the simulation to see whether the lag information contributes anything to the estimation. In total we thus had 10 methods, FOBI, JADE, SOBI, gFOBI and gJADE and the tensorial version of each. Furthermore, for the SOBI-based methods we used the lag set $\mathcal{T}=\{1,2,\ldots,12\}$ and for methods based on gFOBI and gJADE the lag set $\mathcal{T}=\{0,1,\ldots,12\}$, see Matilainen et al. [2015] for previous use of these particular lag sets. For the reasons discussed in the introductory section we expect TSOBI to outperform the other methods in the ARMA case and TgJADE to do the same in the stochastic volatility case.

To make the results of the vector-valued and tensor-valued methods comparable Kronecker products, $\hat{\boldsymbol{\Gamma}}{}^{3}\otimes\hat{\boldsymbol{\Gamma}}{}^{2}\otimes\hat{\boldsymbol{\Gamma}}{}^{1}$, of the 3 unmixing matrix functionals produced by each of the tensor methods were taken. Thus both the Kronecker products and the matrices estimated by the vector methods estimate the inverse of the same $12\times 12$ matrix $\boldsymbol{A}_{3}\otimes\boldsymbol{A}_{2}\otimes\boldsymbol{A}_{1}$. We then measured how close each estimate is to the true inverse via the minimum distance index [Ilmonen et al., 2010]:

[TABLE]

where $p$ is now 12 and $\hat{\boldsymbol{\Gamma}}$ is the estimated Kronecker product or unmixing matrix. The value of the index can be shown to vary from 0 to 1, the former indicating perfect separation.

The mean minimum distance indices (MDI) over different combinations of mixing, method, setting and series length are depicted in Figure 2. Based on the plots, the overall trend is that the vector methods offer no competition to the tensorial ones, the only exception being SOBI which in the ARMA setting outmatches TJADE (and also TFOBI in the case of normal mixing). Furthermore, in the ARMA setting the methods exploiting temporal dependence surpass the ones assuming the independence of observations and in the stochastic volatility setting the best performance is by a clear margin given by a method utilizing temporal dependence, TgJADE. Otherwise the plots offer no surprises: TSOBI, operating solely on autocovariances, outperforms the other methods in the ARMA setting and the FOBI and JADE-based tensor methods dominate in the stochastic volatility setting where no second-order information is present. In the latter case TJADE is almost as good as TgJADE indicating that most of the separation information is in the stochastic volatility setting contained in the higher marginal moments of the series and not in the temporal dependence. Of interest is also the better performance of gFOBI relative to the other non-SOBI vector methods in the ARMA setting, the same phenomenon that was observed also in Matilainen et al. [2015]. The figure also clearly shows that the tensor methods are not affine equivariant - but the dependence on the mixing matrix at least in this case seems to be negligible, especially when compared to the performance of the affine equivariant vectorized versions of the methods.

5.2 Application: fMRI-data

We next consider an application of the methods to task-based fMRI-data simulated using the R-package neuRosim [Welvaert et al., 2011]. The package offers a variety of options for simulating 3-dimensional, noisy, single-subject fMRI-measurements and as our setting we use a particular subset of the settings used in Virta et al. [2016d]. The simulated tensor-valued time series are of length $T=100$ and the individual 3D-frames are of size $64\times 64\times 64$ containing two activation signals residing in disjoint regions, the estimation of which is our objective. The activation signals shown in Figure 3 are convolutions of a stimulus function and a haemodynamic response function (HRF), with two different choices for the latter, gamma and double-gamma. See the section IV.A of Virta et al. [2016d] for more elaborate description of the simulation settings, the simulation code and images of the activation regions.

Our objective is then to estimate the activation signals in Figure 3 and for this we compared a total of 8 methods: TFOBI, TJADE, TSOBI, TgFOBI, TgJADE and the vector-valued versions of all non-JADE-based methods, for which the data cubes were first vectorized. The reason for excluding JADE and gJADE is that their computation was not feasible on a computer with 24 cores and 64GB memory in reasonable time. As is standard in the processing of fMRI-data, prior to applying the methods we first reduced the dimensions of the observed series using principal component-based methods. For the tensor-valued methods tensorial PCA [Virta et al., 2016d] was used to reduce to the case of $6\times 6\times 2$ tensors and for the vector-valued methods SVD was used to compress the observations into vectors of length 72, which are then of comparable dimensionality with the reduced tensors.

We did 1000 replications under both described HRF-settings and for each replication recorded the highest absolute correlations with each of the true signals found amongst the obtained components. Thus the value of one indicates that the method succeeded exactly in extracting the corresponding signal. The resulting boxplots are shown in Figure 4, with the y-axis obeying arctanh-scale (Fisher’s z-transformation) to better emphasize the upper part of its range. The results clearly indicate that apart from SOBI the vectorial methods offer no competition to the tensorial methods. And while the first activation signal is for both HRFs estimated almost equally well by all tensor methods their performances differ in finding the second signal, TSOBI triumphing over the other methods, which is not too surprising as the signals of interest have rather regular patterns and do not really exhibit features of stochastic volatility. Note that the results for TFOBI, TJADE and FOBI in Figure 4 match those obtained in the same settings in Virta et al. [2016d].

5.3 Application: video processing

As described in the introduction, one of the more natural examples of tensor-valued time series is a video clip. A color video with resolution $h\times w$, 3 color channels and $T$ frames can be thought of as the sequence of realizations of a tensor-valued time series $\mathbb{X}_{t}\in\mathbb{R}^{h\times w\times 3}$ for $t=1,\ldots,T$. As our third example we consider the WaterSurface video available online at http://pages.cs.wisc.edu/~jiaxu/projects/gosus/supplement that has been used for background subtraction in Toyama et al. [1999], Li et al. [2004]. The same video was also used for outlier detection in Rousseeuw et al. [2016] and they have the preprocessed video available at https://wis.kuleuven.be/stat/robust/software. The video data tensor consists of a total of $T=633$ color frames of size $128\times 160$ depicting a beach, sea and a tree, with a man entering the scene from left and passing the tree during frames 480-500 and staying in the picture for the rest of the clip. As an illustration Figure 5 shows frames 1 and 520 representing the typical view of the video without and with the man.

Our objective with the video is similar as in Rousseeuw et al. [2016], trying to detect the time window during which the man is in the scene. Assuming that the model (9) holds in this case, the relative staticity of the video before the man’s entrance makes it reasonable to assume that the change point is also visible in the extracted source signals in $\mathbb{Z}_{t}$. For demonstration purposes we considered for this example only TSOBI and extracted the four components with the highest kurtoses and the four components with the lowest kurtoses shown in Figure 6. We chose kurtosis as a criterion for selecting the most interesting components as bimodal distributions are often characterized by extremal kurtosis values. The change point is now visible in the first three source series on the left-hand side of Figure 6, the latter two series having their peaks only after the man has already passed the tree. The fourth series on the left corresponds to a moment when the man has reached the edge of the screen and starts to slowly turn around. Note that contrary to the outlier detection method in Rousseeuw et al. [2016] which flagged all frames after the change point as outlying, TSOBI produced signals that identify the frames during which big changes occur. Interestingly, the first two source series on the right also capture some systematic component in the video although exactly which is not evident from the video. Finally, the running time of the TSOBI-algorithm for the considered video tensor of size $128\times 160\times 3\times 633$ is roughly ten minutes, making it a suitable method also for large-scale problems. Note that using vector-valued BSS methods without resorting to some form of a priori dimension reduction is not possible in this case, the vectorized data set having far more variables ($p=61440$) than observations ($T=633)$.

6 Discussion

Higher order structures are nowadays an increasingly popular form of data and manifest, for example, in the form of video data. One way to treat such observations is as a time series where at each time point we observe a tensor of the same size which, although seemingly natural, as a viewpoint is still largely unexplored in the literature. However, with the increasingly high-dimensional structures we also face the problem of separating the actual information content from the noise. In multivariate time series analysis a standard solution is blind source separation, a technique which in this paper was extended to tensor-valued time series.

Our approach is model-based and we first defined the tensor blind source separation model, a natural extension of the corresponding vector-valued model, where the objective is to recover a latent tensor-valued time series subjected to multimodal linear transformation. Next, we presented tensorial versions of three existing blind source separation methods, SOBI, gFOBI and gJADE, and proved that each of them is capable of estimating the latent source series under specific assumptions. The methods operate respectively on second-order, fourth-order and joint lagged fourth-order moments making the SOBI-based method specifically suited for the separation of linear processes and the last two methods suited to cases where no second-order information exists.

Simulation studies showed that under the model the proposed methods were highly superior not only to the often used combination of vectorizing the tensors and using standard multivariate methods but also to two tensorial dimension reduction methods designed for i.i.d data, TFOBI and TJADE. Hence, if the data exhibits a natural tensor structure and temporal correlation both should also be taken into account in the analysis and processing of the data. In addition, two experiments conducted respectively on simulated fMRI-data and real video data revealed that the proposed methods are also adept in condensing information and extracting signals of importance from high-dimensional tensor-valued time series.

To also theoretically evaluate the gain in performance, our future plans include the investigation of whether asymptotical results like those obtained for SOBI in Miettinen et al. [2016] can be obtained also for the proposed tensor methods, especially TSOBI. Furthermore, we also plan to extend standard time series BSS methods that allow relaxing the stationarity assumption, such as Choi and Cichocki [2000b, a], Nordhausen [2014], to the tensorial case. Another interesting problem is the choice of the most important source components in applications. In Section 5 we used absolute correlation and extremal kurtosis as criteria but several others, such as the number of change points, could be used as well. Also this we plan to explore further in our future work.

7 Acknowledgements

The authors would like to thank the two anonymous referees for their helpful comments and insightful suggestions.

Appendix A The proofs and simulation distributions

The proof of Theorem 2.

Note first that equation (17) in Virta et al. [2016b] says that the $m$-flattening $(\boldsymbol{X}^{st}_{t})^{(m)}$ is equal to $d\cdot\boldsymbol{U}_{m}\boldsymbol{Z}^{(m)}_{t}(\boldsymbol{U}_{m+1}\otimes\cdots\otimes\boldsymbol{U}_{r}\otimes\boldsymbol{U}_{1}\otimes\cdots\otimes\boldsymbol{U}_{m-1})$, where $d$ is the constant of proportionality and the Kronecker product as a product of orthogonal matrices is itself an orthogonal matrix.

Theorem 1 then implies that $\mathbb{X}^{st}_{t}\propto\mathbb{Z}_{t}\odot_{1}\boldsymbol{U}_{1}\cdots\odot_{r}\boldsymbol{U}_{r}$ which, when plugged in into $\boldsymbol{\Sigma}_{\tau}^{m}(\mathbb{X}^{st}_{t})$ and $\boldsymbol{B}_{\tau}^{m}(\mathbb{X}^{st}_{t})$ in combination with the previous result, (T4) and (T5) directly yields the result for the parts i) and ii).

To prove part iii), it is straightforward to show that for any lags $\tau_{1},\tau_{2},\tau_{3},\tau_{4}$ and indices $i,j=1,\ldots,p_{m}$ we have

[TABLE]

where $(\boldsymbol{U}_{m})_{kl}=u^{(m)}_{kl}$. By the assumption we have for any fixed $\tau,i,j$ that $\boldsymbol{C}^{m}_{\tau ij}(\mathbb{Z}_{t})=\boldsymbol{D}^{m}_{\tau ij}$ is diagonal and the matrix $\boldsymbol{C}_{\tau ij}^{m}(\mathbb{X}^{st}_{t})$ can thus be written as

[TABLE]

where $\boldsymbol{\Sigma}^{m}_{0}(\mathbb{Z}_{t})$ is diagonal and by element-wise inspection the quantity inside the square brackets can be showed to be diagonal as well. Thus the right-hand side expression is the eigendecomposition of $\boldsymbol{C}_{\tau ij}^{m}(\mathbb{X}^{st}_{t})$ proving the claim. ∎

The models used in the simulation study.

In the ARMA setting the 12 used time series models were: AR(1) with $\phi=0.9$; AR(1) with $\phi=-0.9$; MA(2) with $\boldsymbol{\theta}=(0.5,-0.5)$; AR(2) with $\boldsymbol{\phi}=(-0.5,-0.3)$; ARMA(4,2) with $\boldsymbol{\phi}=(0.5,-0.3,0.1,-0.1)$ and $\boldsymbol{\theta}=(0.7,-0.3)$; ARMA(2,4) with $\boldsymbol{\phi}=(-0.7,0.1)$ and $\boldsymbol{\theta}=(0.9,0.3,0.1,-0.1)$ and MA(5), MA(10), MA(20), MA(30), MA(40) and MA(50) with parameters sampled independently from the uniform distribution on $(-1,1)$.

In the SV setting the used models consisted of both GARCH and stochastic volatility models, SV($\mu,\phi,\sigma,\nu$), where $\mu$, $\phi$ and $\sigma$ are respectively the mean, AR-parameter and the volatility of the latent AR(1)-process and $\nu$ is the number of degrees of freedom of the innovation $t$-distribution, see Kastner and Frühwirth-Schnatter [2014]. The six used stochastic volatility models were then SV($-10,0.98,0.2,\infty$), SV($-5,-0.98,0.2,10$), SV($-10,0.7,0.7,\infty$), SV($-5,-0.70,0.7,10$), SV($-9,0.20,0.01,\infty$) and SV($-9,-0.20,0.01,10$). In addition, the six used GARCH models were: ARCH(1) with $\alpha=0.7$; GARCH(1,1) with $\alpha=0.2$ and $\beta=0.2$; GARCH(1,1) with $\alpha=0.1$ and $\beta=0.8$; ARCH(4) with $\boldsymbol{\alpha}=(0.20,0.10,0.05,0.01)$; GARCH(3,1) with $\boldsymbol{\alpha}=(0.05,0.03,0.01)$ and $\beta=0.5$ and ARCH(10) with $\boldsymbol{\alpha}=(0.20,0.14,0.12,0.10,0.05,0.05,0.04,0.03,0.02,0.01)$, see Bollerslev [1986] for more details on ARCH and GARCH models. ∎