Streaming data recovery via Bayesian tensor train decomposition

TL;DR

This paper introduces a Bayesian tensor train decomposition method for recovering high-order streaming data, effectively handling incomplete and noisy data by leveraging a streaming variational Bayes approach.

Contribution

It integrates the tensor train format into Bayesian tensor decomposition for streaming data, enabling efficient recovery of high-order, incomplete, and noisy data streams.

Findings

Outperforms existing Bayesian tensor methods in accuracy

Effective with high-order, incomplete, and noisy streaming data

Validated on synthetic and real-world datasets

Abstract

In this paper, we study a Bayesian tensor train (TT) decomposition method to recover streaming data by approximating the latent structure in high-order streaming data. Drawing on the streaming variational Bayes method, we introduce the TT format into Bayesian tensor decomposition methods for streaming data, and formulate posteriors of TT cores. Thanks to the Bayesian framework of the TT format, the proposed algorithm (SPTT) excels in recovering streaming data with high-order, incomplete, and noisy properties. The experiments in synthetic and real-world datasets show the accuracy of our method compared to state-of-the-art Bayesian tensor decomposition methods for streaming data.

Equations

$$x_{\mathbf{j}}\approx a_{\mathbf{j}}=\prod _{d=1}^D{\mathcal{G}}_{j_d}^{(d)},\mathcal{X}\approx \mathcal{A}=⟨⟨{\mathcal{G}}^{(1)},\dots ,{\mathcal{G}}^{(D)}⟩⟩$$

$$\mathcal{L}(\mathcal{A})=\Vert {\mathcal{X}}_{\Omega }-{\mathcal{A}}_{\Omega }{\Vert }_F^2$$

$$⟨\mathcal{X},\mathcal{Y}⟩:=\sum _{\mathbf{j}}{x}_{\mathbf{j}}{y}_{\mathbf{j}}.$$

$$p(\{\boldsymbol{\mathscr{G}}^{(d)}\}_{d=1}^{D})=\prod_{d=1}^{D}\prod_{k=1}^{r_{d-1}}\prod_{l=1}^{r_{d}}\mathcal{N}\left(\boldsymbol{\mathscr{G}}^{(d)}_{k,:,l}\bigg{|}\mathbf{m}_{\boldsymbol{\mathscr{G}}_{k,:,l}^{(d)}},v\mathbf{I}\right),$$

$$p(x_{\mathbf{j}}|\{\boldsymbol{\mathscr{G}}^{(d)}\}_{d=1}^{D},\tau)=\mathcal{N}\left(x_{\mathbf{j}}\big{|}\prod_{d=1}^{D}\boldsymbol{\mathscr{G}}^{(d)}_{j_{d}},\tau^{-1}\right),$$

$$p(\tau |{\alpha }_0,{\beta }_0)=\mathrm{Gam}(\tau |{\alpha }_0,{\beta }_0).$$

$$p(x|\alpha ,\beta )=\frac{x^{\alpha -1}{e}^{-\beta x}{\beta }^{\alpha }}{\Gamma (\alpha )},$$

$$p\left(x_{\mathbf{j}\in \Omega },\{{\mathcal{G}}^{(d)}{\}}_{d=1}^D,\tau \right)$$

$$\displaystyle\cdot\prod_{d=1}^{D}\prod_{k=1}^{r_{d-1}}\prod_{l=1}^{r_{d}}\left(\boldsymbol{\mathscr{G}}^{(d)}_{k,:,l}\bigg{|}\mathbf{m}_{\boldsymbol{\mathscr{G}}_{k,:,l}^{(d)}},v\mathbf{I}\right)$$

$$\displaystyle\cdot\prod_{\mathbf{j}\in\Omega}\mathcal{N}\left(x_{\mathbf{j}}\big{|}\prod_{d=1}^{D}\boldsymbol{\mathscr{G}}^{(d)}_{j_{d}},\tau^{-1}\right).$$

$$p(\theta |B_1)\propto p(\theta )p(B_1|\theta ),$$

$$p(\theta |\tilde{B}\cup B_t)\propto p(\theta |\tilde{B})p(B_t|\theta ),$$

$$q_t(\theta )=q_t(\{{\mathcal{G}}^{(d)}{\}}_{d=1}^D,\tau )=q_t(\tau )\prod _{d=1}^D\prod _{k=1}^{r_{d-1}}\prod _{l=1}^{r_d}\prod _{j_d}^{N_d}{q}_t\left({\mathcal{G}}_{k,j_d,l}^{(d)}\right).$$

$$\mathcal{L}(q_t(\theta ))={\mathbb{E}}_{q_t}\left[\log \left(\frac{q_{t-1}(\theta )p(B_t|\theta )}{q_t(\theta )}\right)\right].$$

$$\log q_t({\theta }_i)={\mathbb{E}}_{q_t(\theta \setminus {\theta }_i)}\left[\log q_{t-1}(\theta )p(B_t|\theta )\right]+\mathrm{const},$$

$$\log q_t({\mathcal{G}}_{k,j_d,l}^{(d)})={\mathbb{E}}_{q_t(\theta \setminus {\mathcal{G}}_{k,j_d,l}^{(d)})}\left[\log q_{t-1}\left(\theta \right)p(B_t|\theta )\right]+\mathrm{const}.$$

$$q_t({\mathcal{G}}_{k,j_d,l}^{(d)})=\mathcal{N}({\mathcal{G}}_{k,j_d,l}^{(d)}|{\hat{m}}_{{\mathcal{G}}_{k,j_d,l}^{(d)}},{\hat{v}}_{{\mathcal{G}}_{k,j_d,l}^{(d)}}).$$

$${\hat{v}}_{{\mathcal{G}}_{k,j_d,l}^{(d)}}={\left(v_{{\mathcal{G}}_{k,j_d,l}^{(d)}}^{-1}+\tau \sum _{\mathbf{j}\in {\Omega }_t^{j_d}}\mathbb{E}\left[b_{(k-1)r_{d-1}+k}^{(<d{)}^T}\right]\mathbb{E}\left[b_{(l-1)r_d+l}^{(>d)}\right]\right)}^{-1},$$

$${\hat{m}}_{{\mathcal{G}}_{k,j_d,l}^{(d)}}$$

$$\sum _{k^{\prime }=1,k^{\prime }\ne k}^{r_{d-1}}\sum _{l^{\prime }=1,l^{\prime }\ne l}^{r_d}\mathbb{E}\left[b_{(k-1)r_{d-1}+k^{\prime }}^{(<d{)}^T}\right]m_{{\mathcal{G}}_{k^{\prime },j_d,l^{\prime }}^{(d)}}\mathbb{E}\left[b_{(l-1)r_d+l^{\prime }}^{(>d)}\right])),$$

$$\mathbb{E}\left[e_k^{(<d{)}^T}\right]=\prod _{i=1}^{d-1}\mathbb{E}\left[{\mathcal{G}}_{j_i}^{(i)}\right],\mathbb{E}\left[e_k^{(>d)}\right]=\prod _{i=d}^D\mathbb{E}\left[{\mathcal{G}}_{j_i}^{(i)}\right],$$

$$\mathbb{E}\left[b^{(<d{)}^T}\right]=\prod _{i=1}^{d-1}\mathbb{E}\left[{\mathcal{G}}_{j_i}^{(i)}\otimes {\mathcal{G}}_{j_i}^{(i)}\right],\mathbb{E}\left[b^{(>d)}\right]=\prod _{i=d+1}^D\mathbb{E}\left[{\mathcal{G}}_{j_i}^{(i)}\otimes {\mathcal{G}}_{j_i}^{(i)}\right].$$

$$\mathbb{E}\left[{\mathcal{G}}_{j_i}^{(i)}\otimes {\mathcal{G}}_{j_i}^{(i)}\right]=\mathbb{E}\left[{\mathcal{G}}_{j_i}^{(i)}\right]\otimes \mathbb{E}\left[{\mathcal{G}}_{j_i}^{(i)}\right]+\mathbb{E}\left[\left({\mathcal{G}}_{j_i}^{(i)}-\mathbb{E}\left[{\mathcal{G}}_{j_i}^{(i)}\right]\right)\otimes \left({\mathcal{G}}_{j_i}^{(i)}-\mathbb{E}\left[{\mathcal{G}}_{j_i}^{(i)}\right]\right)\right],$$

$$\mathbb{E}\left[{\mathcal{G}}_{k,j_i,l}^{(i)}{\mathcal{G}}_{k^{\prime },j_i,l^{\prime }}^{(i)}\right]=\mathbb{E}\left[{\mathcal{G}}_{k,j_i,l}^{(i)}\right]\mathbb{E}\left[{\mathcal{G}}_{k^{\prime },j_i,l^{\prime }}^{(i)}\right]+v_{{\mathcal{G}}_{k,j_i,l}^{(i)}}\delta \left(k-k^{\prime }\right)\delta \left(l-l^{\prime }\right).$$

$$\mathbb{E}\left[{\mathcal{G}}_{k,j_i,l}^{(i)}{\mathcal{G}}_{k^{\prime },j_i,l^{\prime }}^{(i)}\right]=\mathbb{E}\left[{\mathcal{G}}_{k,j_i,l}^{(i)}\right]\mathbb{E}\left[{\mathcal{G}}_{k^{\prime },j_i,l^{\prime }}^{(i)}\right]+v_{{\mathcal{G}}_{k,j_i,l}^{(i)}}\delta \left(k-k^{\prime }\right)\delta \left(l-l^{\prime }\right),$$

$$\mathbb{E}{\left[{\mathcal{G}}_{j_i}^{(i)}\otimes {\mathcal{G}}_{j_i}^{(i)}\right]}_{p,q}={\left(\mathbb{E}\left[{\mathcal{G}}_{j_i}^{(i)}\right]\otimes \mathbb{E}\left[{\mathcal{G}}_{j_i}^{(i)}\right]\right)}_{p,q}+{\mathcal{V}}_{p,j_i,q}^{(i)},$$

$${\mathcal{V}}_{j_d}^{(i)}:=\mathbb{E}\left[\left({\mathcal{G}}_{j_i}^{(i)}-\mathbb{E}\left[{\mathcal{G}}_{j_i}^{(i)}\right]\right)\otimes \left({\mathcal{G}}_{j_i}^{(i)}-\mathbb{E}\left[{\mathcal{G}}_{j_i}^{(i)}\right]\right)\right].$$

$$\mathbb{E}\left[{\mathcal{G}}_{j_i}^{(i)}\otimes {\mathcal{G}}_{j_i}^{(i)}\right]=\mathbb{E}\left[{\mathcal{G}}_{j_i}^{(i)}\right]\otimes \mathbb{E}\left[{\mathcal{G}}_{j_i}^{(i)}\right]+{\mathcal{V}}_{j_i}^{(i)}.$$

$$\log q_t(\tau )={\mathbb{E}}_{q_t(\theta \setminus \tau })\left[\log q_{t-1}\left(\theta \right)p(B_t|\theta )\right]+\mathrm{const},$$

$$q_t(\tau )=\mathrm{Gam}(\tau |\hat{\alpha },\hat{\beta }),$$

$$\hat{\alpha }=\alpha +\frac{|{\Omega }_t|}{2},$$

$$\hat{\beta }=\beta +\frac{1}{2}\mathbb{E}\left[\sum _{\mathbf{j}\in {\Omega }_t}{\left(x_{\mathbf{j}}-\prod _{d=1}^D{\mathcal{G}}_{j_d}^{(d)}\right)}^2\right],$$

Full text

††journal: Applied Mathematics and Computation

1 Introduction

Multiway data or multidimensional arrays, represented by tensors, are ubiquitous in real-world applications, such as recommendation systems [1, 2], computer version [3], medical imaging [4] and uncertainty quantification [5]. How to process, analyze and utilize such high-volume tensor data is a fundamental problem in machine learning [6, 7], data mining [8] and signal processing [9].

Effective numerical techniques, such as CANDECOMP/PARAFAC (CP) decomposition [10, 11, 12] and Tucker decomposition [13, 14] are the most commonly used tensor decomposition approaches and have been proposed to compress full tensors and to obtain their low-rank representations. CP decomposition approximates a tensor by a sum of rank-one tensors, while Tucker decomposition decomposes a tensor into a core tensor and several factor matrices. Since CP decomposition can be seen as a special case of Tucker decomposition [15], Tucker decomposition is more flexible than CP decomposition. However, due to the existence of a core tensor, Tucker decomposition also brings challenges in both modeling and computation. In this paper, we mainly focus on Tensor Train (TT) decomposition [16], which combines the advantages of CP and Tucker decomposition, because it provides a space-saving model called TT format while preserving the representation power.

This paper is interested in the decomposition of streaming data. Due to the stress on database capacity and privacy, streaming data is generated continuously by different of data sources and in small sizes, such as log files from web application [17] and information from social networks [18]. Recently, several works decompose fast streaming data, e.g. POST [19] and BASS-Tucker [20]. POST and BASS-Tucker are probabilistic streaming tensor decomposition algorithm based on Bayesian formulation of CP decomposition and Tucker decomposition respectively. Although both of POST and BASS-Tucker are powerful streaming decomposition approach and can adapt to a large amount of applications, the latent factors (CP factors) in POST are lack of representation power and BASS-Tucker is ambitious to execute and comprehensive due to the complex reducing model and high time complexity.

To address these challenges, we develop SPTT, a Streaming Probabilistic TT decomposition approach. The main contributions of this work are three-folds: First, we propose a Bayesian tensor decomposition method based on TT decomposition, which can be adopted to analysis high-order incomplete noisy streaming data; Second, motivated by the Bayesian TT representation for batch-mode data in [21], we introduce a probabilistic TT model and exploit the streaming variational inference (SVI) method on it for streaming data analysis; Third, both synthetic and real-world data are used in demonstrating the accuracy of our algorithm.

The paper is organized as follows. We review the tensor basics and tensor train decomposition in section 2, and in 3 we present our SPTT algorithm with streaming variational inference. In section 4, we demonstrate the accuracy of our algorithm in both synthetic and real-worlds dataset. Finally, in section section 5, we make some conclusions.

2 Background

2.1 Tensor basics and notation

In this paper, we denote scalar, vector, matrix and tensor variables, respectively, by lowercase letters (e.g. $x$), boldface lowercase letters (e.g. $\mathbf{x}$), boldface capital letters (e.g. $\mathbf{X}$), and boldface Euler script letters (e.g. $\boldsymbol{\mathscr{X}}$). For a $D$-th order tensor $\boldsymbol{\mathscr{X}}\in\mathbb{R}^{N_{1}\times\cdots\times N_{D}}$, we denote the $\mathbf{j}:=(j_{1},\cdots,j_{D})$-th element of $\boldsymbol{\mathscr{X}}$ by $x_{\mathbf{j}}$, where $j_{d}=1,\ldots,N_{d}$ for $d=1,\ldots,D$.

2.2 Tensor train decomposition

Given TT-ranks $r_{k}$ in vector $\mathbf{r}=(r_{1},\ldots,r_{D})$ with $r_{1}=r_{D}=1$, each element of a $D$-th order tensor $\boldsymbol{\mathscr{X}}$ is tensorized by TT decomposition

[TABLE]

where $\boldsymbol{\mathscr{G}}^{(d)}_{j_{d}}\in\mathbb{R}^{r_{d-1}\times r_{d}}$ is the $j_{d}$-th slice of TT-cores $\boldsymbol{\mathscr{G}}^{(d)}\in\mathbb{R}^{r_{d-1}\times N_{d}\times r_{d}}$, for $j_{d}=1,\ldots,N_{d}$, $d=1,\ldots,D$, and $a_{\mathbf{j}}$ is the $\mathbf{j}$-th element of tensor $\boldsymbol{\mathscr{A}}$ and approximates $x_{\mathbf{j}}$. In other words, TT decomposition is to decompose a $D$-th order tensor $\boldsymbol{\mathscr{X}}$ into a sequence of three-way factor tensors, collected by a TT-format tensor $\boldsymbol{\mathscr{A}}$. Figure 1 illustrate the TT decomposition.

To find the optimal TT-cores under partial observations, TT decomposition minimizes a mean square loss,

[TABLE]

where $\boldsymbol{\mathscr{X}}_{\Omega}$ and $\boldsymbol{\mathscr{A}}_{\Omega}$ are vectors containing $x_{\mathbf{j}}$ and $a_{\mathbf{j}}$, respectively, $\mathbf{j}\in\Omega$. Here $\Omega$ is the index set of observed tensor entries. We can use alternating least squares [22] to update each TT core alternatively, given all the other fixed. In addition, gradient-based approaches [23] can also be used to minimize the loss function.

Instead of explicitly storing the original tensor, one can choose to maintain and use the TT-cores because of the convenience and effectiveness of TT-cores. Assuming that $r_{d}=R,\,d=2,\dots.D-1$, and that $N_{d}=N,\,d=1,\dots,D$, the storage of the tensor $\boldsymbol{\mathscr{X}}$ in the standard form is $O(N^{D})$, while the storage in TT-format $\boldsymbol{\mathscr{A}}$ is $O(NDR^{2})$, which is linear to the tensor order. Consider two tensors $\boldsymbol{\mathscr{X}},\boldsymbol{\mathscr{Y}}\in\mathbb{R}^{N_{1}\times\cdots\times N_{D}}$, the inner product is defined as

[TABLE]

Similarly, using the TT-format of $\boldsymbol{\mathscr{X}},\boldsymbol{\mathscr{Y}}$, one can implement the inner product with $O(DNR^{3})$ operations rather than $O(N^{D})$ (Cf. Algorithm 4 in [16]). Using the inner product, the Frobenius norm $\|\boldsymbol{\mathscr{X}}\|_{F}=\sqrt{\langle\boldsymbol{\mathscr{X}},\boldsymbol{\mathscr{X}}\rangle}$ and the distance between two tensors $\|\boldsymbol{\mathscr{X}}-\boldsymbol{\mathscr{Y}}\|_{F}$ can be efficiently calculated. Besides, when TT-ranks are small enough, the TT-format is also efficient for other basic linear algebra operations, such as matrix-by-vector products, scalar multiplications, see [16] for details.

3 Bayesian model for streaming tensor train decomposition

In this section, we first develop the probabilistic model of TT decomposition, where the Gaussian noise is included. After that, streaming data is introduced and a posterior inference for TT-cores is developed by streaming variational inference (SVI) [24]. Finally, we present the streaming probabilistic tensor train decomposition (SPTT) algorithm.

3.1 Probabilistic modeling of tensor train decomposition

The standard TT decomposition, like [22, 23], use the point estimation to approximate the TT-cores and is not capable of evaluating the uncertainty, which can vary among each slice in a TT core. To overcome this problem, inspired by a Bayesian TT representation for batch-mode data in [21], we introduce a Bayesian TT model for streaming data.

We first specify $\{\boldsymbol{\mathscr{G}}^{(d)}\}_{d=1}^{D}$ from a Gaussian prior,

[TABLE]

which indicates that each fiber $\boldsymbol{\mathscr{G}}^{(d)}_{k,:,l}$ admits a Gaussian distribution with mean $\mathbf{m}_{\boldsymbol{\mathscr{G}}_{k,:,l}}^{(d)}$, and covariance matrix $v\mathbf{I}$. $v$ is a scalar and controls the flatness of the Gaussian prior.

Given the TT-cores $\{\boldsymbol{\mathscr{G}}^{(d)}\}_{d=1}^{D}$, the likelihood of an observed entry value $x_{\mathbf{j}}$ is

[TABLE]

where $\tau$ is the inverse of the noise variance. Since the Gamma distribution, because of its non-negative and long tail, provides a good model for $\tau$, we assign a Gamma prior over $\tau$

[TABLE]

where $\alpha_{0}$ and $\beta_{0}$ are hyperparameters of the Gamma distribution. The density function of a Gamma distribution is

[TABLE]

where $\Gamma(\alpha)$ is a Gamma function.

The joint probability is then given by combining (4) with (5)-(6),

[TABLE]

3.2 Bayesian streaming inference

We introduce SPTT, our posterior inference algorithm for streaming data, based on TT decomposition in subsection 3.1 and the streaming variational inference (SVI) [24]. Assume that the observed tensor entries are provided in a series of batches, $\{B_{1},B_{2},\ldots\}$, we aim to conduct the posterior inference for TT-cores $\{\boldsymbol{\mathscr{G}}^{(d)}\}_{d=1}^{D}$ and the inverse of noise variance $\tau$ upon receiving each data batch $B_{t}$, without revisiting the previous batches $\tilde{B}=\{B_{1},\ldots,B_{t-1}\}$.

Assume the prior distribution of $\mathbf{\theta}=\{\{\boldsymbol{\mathscr{G}}^{(d)}\}_{d=1}^{D},\tau\}$ is $p(\mathbf{\theta})$. Denote $\Omega_{t}$ is the index set of $t$-th data batch $B_{t}$. For the first data batch $B_{1}=\{x_{\mathbf{j}}\}_{\mathbf{j}\in\Omega_{1}}$, the posterior $p(\mathbf{\theta}|B_{1})$ can be calculated via Bayes’ rule:

[TABLE]

where $p(B_{1}|\mathbf{\theta})=\prod_{\mathbf{j}\in\Omega_{1}}p(x_{\mathbf{j}}|\mathbf{\theta})$. Similarly, when receiving the newly data batch $B_{t}=\{x_{\mathbf{j}}\}_{\mathbf{j}\in\Omega_{t}}$, the incremental posterior is given by

[TABLE]

The proportional relationship in (9) motivates us to adopt a recursive process of streaming inference for $\mathbf{\theta}$. However, the multilinear product of TT-cores $\{\boldsymbol{\mathscr{G}}^{(d)}\}_{d=1}^{D}$ in (5) results in the posterior update intractable. To overcome the problem, we introduce a factorized variational posterior distribution $q_{t-1}(\mathbf{\theta})=\prod_{i}q_{t-1}(\mathbf{\theta}_{i})$ to approximate $p(\mathbf{\theta}|\tilde{B})$. With receiving new data batch $B_{t}$, we can obtain a new distribution $\tilde{q}_{t}(\mathbf{\theta})=q_{t-1}(\mathbf{\theta})p(B_{t}|\mathbf{\theta})$, which is proportional to $p(\mathbf{\theta},\tilde{B}\cup B_{t})$. We approximate posterior $p(\mathbf{\theta}|\tilde{B}\cup B_{t})$ with $q_{t}(\mathbf{\theta})$,

[TABLE]

based on the only assumption that each entries in the TT-cores and the inverse of noise variance are independent. To be noticed, $q_{t}(\theta)$ is also in a factorized form as $q_{t-1}(\mathbf{\theta})$, and can be initialized with $q_{0}(\mathbf{\theta})=p(\mathbf{\theta})$. From the variational inference framework [25], we minimize the Kullback-Leibler (KL) divergence between $\frac{1}{C}\tilde{q}_{t}(\mathbf{\theta})$ and $q_{t}(\mathbf{\theta})$, i.e. $\operatorname{KL}(\frac{1}{C}\tilde{q}_{t}(\mathbf{\theta}||q_{t}(\mathbf{\theta}))$, where $C$ is a normalization constant. This is equivalent to maximize a variational model evidence lower bound (ELBO) [25],

[TABLE]

Due to the factorized form (10), the particular functional forms of the individual factors $q_{t}(\mathbf{\theta}_{i})$ can be explicitly derived in turn. The update rule for the $i$-th factor based on the maximization of $\mathcal{L}(q_{t}(\mathbf{\theta}))$ is then given by

[TABLE]

where $\mathbb{E}_{q_{t}(\theta\setminus\mathbf{\theta}_{i})}$ denotes an expectation w.r.t the $q_{t}$ distributions over all variables except $\mathbf{\theta}_{i}$.

3.2.1 Posterior distribution of TT-cores

Upon the factorized form (10), the inference of an element of TT-cores $\boldsymbol{\mathscr{G}}_{k,j_{d},l}^{(d)}$ can be performed by receiving the messages from the $t$-th data batch $B_{t}$ and the updated information of the whole TT-cores, where $k=1,\ldots,r_{d-1}$, $l=1,\ldots r_{d}$, $j_{d}=1,\ldots,N_{d}$. By applying (12), the optimized form of $\boldsymbol{\mathscr{G}}_{k,j_{d},l}^{(d)}$ is given as follows,

[TABLE]

Obviously, the approximation posterior keeps the same form as the updated one, i.e., Gaussian distribution,

[TABLE]

The variance and mean, derived from derivation (13) w.r.t $\boldsymbol{\mathscr{G}}_{k,j_{d},l}^{(d)}$, can be updated by

[TABLE]

[TABLE]

with

[TABLE]

[TABLE]

where $\hat{v}_{\boldsymbol{\mathscr{G}}_{k,j_{d},l}^{(d)}}$ and $\hat{m}_{\boldsymbol{\mathscr{G}}_{k,j_{d},l}^{(d)}}$ are the variance and mean of the element $\boldsymbol{\mathscr{G}}_{k,j_{d},l}^{(d)}$ after the $t$-th data batch, while ${v}_{\boldsymbol{\mathscr{G}}_{k,j_{d},l}^{(d)}}$ and ${m}_{\boldsymbol{\mathscr{G}}_{k,j_{d},l}^{(d)}}$ are the updated variance and mean after the $(t-1)$-th data batch. The subscript $\Omega_{t}^{j_{d}}$ represents the indexes of $t$-th batch data with $d$-th index fixing as $j_{d}$. The expression $\mathbb{E}[\cdot]$ denotes the expectation with respect to all random variables involved. From (17) and (18), we can see that, for $d=1,\ldots,D$, $e^{(<d)^{T}}$ and $b^{(<d)^{T}}$ are row vectors, while $e^{(>d)}$ and $b^{(>d)}$ are column vectors. The subscript of them denotes an element in a vector, e.g. $e^{(<d)^{T}}_{k}$ is the $k$-th element in vector $e^{(<d)^{T}}$. To compute $\mathbb{E}\left[\boldsymbol{\mathscr{G}}^{(i)}_{j_{i}}\otimes\boldsymbol{\mathscr{G}}^{(i)}_{j_{i}}\right]$ in (18), we introduce the following lemma.

Lemma 1**.**

The expectation of Kronecker product $\boldsymbol{\mathscr{G}}^{(i)}_{j_{i}}\otimes\boldsymbol{\mathscr{G}}^{(i)}_{j_{i}}$ in (18) can be calculated by

[TABLE]

and its element-wise form is

[TABLE]

Proof.

From equation (10), we can see any two entries in $\boldsymbol{\mathscr{G}}^{(i)}_{j_{i}}$ are independent,

[TABLE]

where $\delta(x)=1$ if $x=0$ and zero otherwise. Let $p=r_{i-1}(k-1)+l,\ q=r_{i}(k^{\prime}-1)+l^{\prime}$, then

[TABLE]

where

[TABLE]

As a Kronecker-form covariance, $\boldsymbol{\mathscr{V}}_{j_{i}}^{(i)}\in\mathbb{R}^{r_{i-1}^{2}\times r_{i}^{2}}$ consists of block matrices with dimension $r_{i-1}\times r_{i}$ and the $(k,l)$-th block matrix contains only one nonzero element $v_{\boldsymbol{\mathscr{G}}_{k,j_{i},l}^{(i)}}$ at the $(k,l)$-th position. Thus,

[TABLE]

∎

3.2.2 Posterior distribution of inverse of noise variance

Following equation (12), the update for the inverse of noise variance $\tau$ can be derived by rearranging equation as follows,

[TABLE]

from which, we can see the posterior of inverse of noise variance is again a Gamma distribution,

[TABLE]

and the posterior hyperparameters can be updated by

[TABLE]

[TABLE]

where

[TABLE]

From (22), we can see that the shape parameter updates with number of the newly receiving data batch $B_{t}$. The summation term in (23) controls the inverse of noise variance $\tau$ through the rate parameter $\hat{\beta}$. It is seen in (23) that when the model dose not fit the streaming data well, there is an increase in $\hat{\beta}$ and thus a decrease in noise variance since $\mathbb{E}[\tau]=\hat{\alpha}/\hat{\beta}$. Computing (24) is challenging, so we present a lemma as follows.

Lemma 2**.**

Given a set of independent random matrices $\{\boldsymbol{\mathscr{G}}_{j_{d}}^{(d)}\}_{d=1}^{D}$, we assume that $\forall d\in\{1,\ldots,D\},\forall k\in\{1\,\ldots,r_{d-1}\},\forall l\in\{1\,\ldots,r_{d}\}$, the vectors $\{\boldsymbol{\mathscr{G}}_{k,:,l}^{(d)}\}$ are independent, then

[TABLE]

Proof.

For a scalar $\prod_{d=1}^{D}\boldsymbol{\mathscr{G}}^{(d)}_{j_{d}}$, using the mix product property of Kronecker product and conclusion of Lemma 1, we have

[TABLE]

where $\boldsymbol{\mathscr{V}}_{j_{d}}^{(d)}$ is defined in (19). ∎

3.2.3 Algorithm

Since variational Bayesian inference is only guaranteed to converge to a local minimum, a good initialization is important. The inverse of noise variance in (6) is initialized as $\mathbb{E}[\tau]=1$ with $q_{0}(\tau)=\operatorname{Gam}(\tau|\alpha_{0},\beta_{0})$ and $\alpha_{0}=\beta_{0}=10^{-3}$. For the TT-cores $\{\boldsymbol{\mathscr{G}}^{(d)}\}_{d=1}^{D}$ in (4), each fiber $\boldsymbol{\mathscr{G}}^{(d)}_{k,:,l}$ is initialized as a vector $\mathbf{m}_{\boldsymbol{\mathscr{G}}^{(d)}_{k,:,l}}$ and its elements are generated from a standard uniform distribution $U(0,1)$. That is,

[TABLE]

The following relative error of TT-cores between two iterations is defined as stopping criterion for data batch $B_{t}$,

[TABLE]

where $\|\cdot\|_{F}$, is the Frobenius norm defined in (3). $\boldsymbol{\mathscr{G}}^{(d)}(m,t)$ denotes the $d$-th tensor core after $m$-th iteration and $t$-th data batch. For each data batch $B_{t}$, we declare our update scheme converged when the relative error $\epsilon_{\mathscr{G}}$ is less than a tolerance $\operatorname{tol}=10^{-3}$. We stop the SPTT algorithm and output the updated mean and variance of each TT-cores if the final data batch $B_{t_{\max}}$ is used. Our algorithm is summarized in Algorithm 1.

The complexity of the proposed algorithm comes from the update the posterior of $\{\boldsymbol{\mathscr{G}}^{(d)}\}_{d=1}^{D}$ and $\tau$. For simplicity, we assume that all TT-ranks are initially set as $L$, the length in each order is $N$, and the batch size is $S$. In each streaming batch, for the updating of each TT core $\boldsymbol{\mathscr{G}}^{(d)}$, it takes $\mathbf{O}(DL^{2})$ to get $e^{(<d)^{T}}$ and $e^{(>d)}$, $\mathbf{O}(DL^{4})$ to get $b^{(<d)^{T}}$ and $b^{(>d)}$, then $\mathbf{O}(SDL^{2})$ to get the mean and $\mathbf{O}(SDL^{4})$ to get the variance of a TT core in (15) and (16) respectively. Furthermore, each $\tau$ requires $\mathbf{O}(SD(L^{4}+L^{2}))$ to compute. Thus, the overall time complexity is $\mathbf{O}(SDL^{4})$. The space complexity is $\mathbf{O}(NDL^{2})$, which is to store the posterior mean and variance of $\{\boldsymbol{\mathscr{G}}^{(d)}\}_{d=1}^{D}$.

4 Experimental results

In this section, we evaluated our streaming probabilistic tensor train decomposition (SPTT) approach on both synthetic data and real-world applications. We compare our proposed SPTT approach with the following tensor decomposition approaches: (1) POST [19], as streaming Bayesian CP decomposition algorithm. (2) BASS-Tucker [20], a streaming Bayesian sparse Tucker decomposition algorithm. (3) CP-ALS [26], a CP decomposition with alternating least square and EM. (4) CP-WOPT [27], which is a static CP decomposition using conjugate gradient descent. (5) Tucker-ALS [28], Tucker decomposition with alternating least square updates. All results of this paper are obtained in MATLAB (2022a) on a computer with an Intel(R) Core(TM) i7-8650U CPU @ 1.90GHz 2.11 GHz .

4.1 Synthetic data

The synthetic tensor data is generated by the following procedure. Through TT decomposition in (1), we create a true tensor $\boldsymbol{\mathscr{A}}\in\mathbb{R}^{20\times 20\times 20\times 20}$ in TT-format with the TT-ranks $(1,3,3,3,1)$ by sampling the entries of TT-cores $\{\boldsymbol{\mathscr{G}}^{(d)}\}_{d=1}^{4}$ from standard uniform distribution $U(0,1)$, and $\boldsymbol{\mathscr{A}}=\langle\langle\boldsymbol{\mathscr{G}}^{(1)},\ldots,\boldsymbol{\mathscr{G}}^{(4)}\rangle\rangle$. The synthetic observed tensor $\boldsymbol{\mathscr{U}}$ is generated as follows,

[TABLE]

where $\boldsymbol{\mathscr{W}}$ is a Gaussian noise tensor with its element $w_{\mathbf{j}}\sim\mathcal{N}(0,\sigma^{2})$ and $\sigma$ is noise level, controlled by signal-to-noise-ratio (SNR),

[TABLE]

To investigate the ability of our streaming decomposition algorithm, we set the ratio of the number of observed tensor entries to the total tensor entries as $15\%$. To be noticed, the index of observed entries are chosen randomly. We divide $10\%$ observed entries into the test set $T$, remaining $90\%$ of observed entries into the training set $B$. Then we randomly partition the training set $B$ into a stream of small data bathes $\{B_{1},B_{2},\ldots,B_{t_{\max}}\}$. Denote batch size as $S$, then $t_{\max}=|B|/S$, where $|B|$ is the total number of training set. Set the maximum iteration number in Algorithm 1 as $M=100$, the initial TT-cores with TT-ranks as $(1,R,R,R,1)$. For maximum iteration number and initialization of all the completing approaches, we use their default settings.

Denote $\hat{\boldsymbol{\mathscr{A}}}$ as the TT-format tensor formulated by the updated posterior mean of TT-cores. To evaluate the performance of the above tensor decomposition approaches, we use the following relative error $\epsilon$ as reconstruction criterion,

[TABLE]

Next, we compare the relative reconstruction error under different batch size, different initial TT-ranks, i.e. different $R$ and different SNR. We repeat our algorithm $5$ times and calculate the average relative error as $\frac{1}{5}\sum_{i=1}^{5}\epsilon_{i}$ and standard deviation are showed in Fig. 2. In Fig. 2 (a), we fix the rank $R=3$, $\operatorname{SNR}=20$, and examine the performance with different choices of batch sizes $B=\{2^{8},2^{9},2^{10},2^{11}\}$. In Fig. 2 (b), we fix the streaming batch size $S=512$, $\operatorname{SNR}=20$, and show how the predictive performance of each method varied with different rank $R=\{3,4,5,6\}$. In Fig. 2 (c), we fix $S=512$, $R=3$, and see the performance with different noise level $\operatorname{SNR}=\{15,20,25,30\}$. As we can see, SPTT outperforms BASS-Tucker, POST and static decomposition methods in all the cases.

In addition, we evaluate the running predictive performance of SPTT. We fix the batch size $S=512$, $\operatorname{SNR}=20$, and gradually feed the training set $\{B_{1},B_{2},\ldots,B_{t_{\max}}\}$ to SPTT, BASS-Tucker and POST. Two cases ($R=3,\,5$) are considered to test the prediction accuracy after each data batch is used. The corresponding running relative error is shown in Fig. 3. As we can see, our SPTT algorithm beats POST and BASS-Tucker since the first data batch, which demonstrates the accuracy and efficiency of our algorithm.

4.2 Real-world applications

We evaluate the SPPT algorithm on the following three datasets: (1) Alog [29], a three-mode (user, action, resource) tensor, of size $200\times 100\times 200$, including $0.66\%$ observed entries. It represents three-way management operations. (2) Carphone [30], a three-mode (length, width, frames) of video data size $144\times 176\times 180$, contains $0.5\%$ nonzeros entries. (3) Flow injection111Retrieved from www.models.kvl.dk, a three-mode (samples, wavelengths, times) tensor of size $12\times 100\times 89$, is exacted from a flow injection analysis system and includes $21360$ nonzero entries.

Set $M=100$ in Algorithm 1 and the initial TT-cores with TT-ranks as $(1,R,R,R,1)$. For maximum iteration number and initialization of all the completing approaches, we use their default settings. We compare the relative reconstruction error $\epsilon$ defined in (29) under different batch size $S$ and different initial rank $R$. We randomly split the observed entries of Alog into $75\%$ for training and $25\%$ for test and the other two datasets $10\%$ entries for test. For POST, SPTT and BASS-Tucker, we randomly partitioned the training entries into a stream of small batches according to batch size $S$. For static decomposition algorithm CP-ALS, CP-WOPT and Tucker-ALS, we input the whole data. The average relative error and standard deviation are reported in Fig. 4. In Fig. 4 (a)-(c), we fix the batch size $S=512$ and show how the predictive performance of each method varied with different rank $R=\{3,4,5,6\}$. The bigger the rank, the more expensive for SPTT to factorize each batch. It can be seen form Fig. 4 (c), as the initial rank $R$ increase from $3$ to $5$, the relative error of SPTT gradually increase and finally close to the relative error of POST. Hence this setting examined the trade-off between the accuracy and computational complexity. To test the capacity of decomposition algorithm on different batch size, we fix the rank $R=3$, and examine the performance with different choices of batch sizes $B=\{2^{8},2^{9},2^{10},2^{11}\}$ in Fig. 4 (d)-(f). As we can see, SPTT outperforms BASS-Tucker, POST and static decomposition methods in three datasets, especially in dataset Carphone.

We then evaluate the running predictive performance of SPTT. We fix the batch size $S=512$, and continuously feed the training set to SPTT, BASS-Tucker and POST. For the three datasets, we consider two cases, i.e. $R=3$ and $R=5$ to test the prediction accuracy after each data batch is used. The running relative error is shown in Fig. 3. In all the cases, when the number of batches is small, POST or BASS-Tucker may performs better than SPTT. However, as the data streams, SPTT beats POST and BASS-Tucker gradually.

5 Conclusions

In this paper, we propose a novel and efficient algorithm, namely SPTT, to discover latent factors from streaming data from high-order tensor and predict missing values. We introduce a simple probabilistic TT model for data denosing and decomposition. Under the SVI framework, we develop posterior inference for TT-cores and the inverse of noise variance. We demonstrate that our SPTT algorithm can produce accurate predictions on streaming data. For future work, we are interested in develop a sparse strategy based on SPTT to further prevent overfitting and enhance scalability.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (No. 12071291), the Science and Technology Commission of Shanghai Municipality (No. 20JC1414300) and the Natural Science Foundation of Shanghai (No. 20ZR1436200).