Bayesian CP Factorization of Incomplete Tensors with Automatic Rank Determination

Introduction

Tensors (i.e., multiway arrays) provide an effective and faithful representation of the structural properties of data, in particular, when multidimensional data or a data ensemble affected by multiple factors are involved. For instance, a video sequence can be represented by a third-order tensor with dimensionality of $height\times width\times time$; an image ensemble measured under multiple conditions can be represented by a higher order tensor with dimensionality of $pixel\times person\times pose\times illumination$. Tensor factorization enables us to explicitly take into account the structure information by effectively capturing the multilinear interactions among multiple latent factors. Therefore, its theory and algorithms have been an active area of study during the past decade (see e.g., ), and have been successfully applied to various application fields, such as face recognition, social network analysis, image and video completion, and brain signal processing. The two most popular tensor factorization frameworks are Tucker and CANDECOMP/PARAFAC (CP), also known as canonical polyadic decomposition (CPD) .

Most existing tensor factorization methods assume that the tensor is complete, whereas the problem of missing data can arise in a variety of real-world applications. This issue has attracted a great deal of research interest in tensor completion in recent years. The objective of tensor factorization of incomplete data is to capture the underlying multilinear factors from only partially observed entries, which can in turn predict the missing entries. In , CP factorization with missing data was formulated as a weighted least squares problem, termed CP weighted optimization (CPWOPT). A structured CPD using nonlinear least squares (CPNLS) was proposed in . In , geometric nonlinear conjugate gradient (geomCG) based on Riemannian optimization on the manifold of tensors were presented. However, as the number of missing entries increases, tensor factorization schemes tend to overfit the model because of an incorrectly specified tensor rank, resulting in severe deterioration of their predictive performance. In contrast, another popular technique is to exploit a low-rank assumption for recovering the missing entries; the rank minimization can be formulated as a convex optimization problem on a nuclear norm. This technique has been extended to higher order tensors by defining the nuclear norm of a tensor, yielding the tensor completion . Some variants were also proposed, such as a framework based on convex optimization and spectral regularization, which uses an inexact splitting method , and fast composite splitting algorithms (FCSA) . To improve the efficiency, the Douglas-Rachford splitting technique , nonlinear Gauss-Seidal method were also investigated. Recently, the nuclear norm based optimization was also applied to a supervised tensor dimensionality reduction method . An alternative method for tensor completion is to employ adaptive sampling schemes . However, the nuclear norm of a tensor is defined straightforwardly by a weighted sum of the nuclear norm of mode-$n$ matricizations, which is related to multilinear rank rather than CP rank. In addition, these completion-based methods cannot explicitly capture the underlying factors. Hence, a simultaneous tensor decomposition and completion (STDC) method was introduced in which a rank minimization technique was combined with Tucker decomposition . To improve completion accuracy, auxiliary information was also exploited in , which strongly depends on the specific application. It is also noteworthy that the rank minimization based on convex optimization of the nuclear norm is affected by tuning parameters, which may tend to over- or under-estimate the true tensor rank.

It is important to emphasize that our knowledge about the properties of CP rank, defined by the minimum number of rank-one terms in CP decomposition, is surprisingly limited. There is no straightforward algorithm to compute the rank even for a given specific tensor, and the problem has been shown to be NP-complete . The lower and upper bound of tensor rank was studied in . The ill-posedness of the best low-rank approximation of a tensor was investigated in . In fact, determining or even bounding the rank of an arbitrary tensor is quite difficult in contrast to the matrix rank , and this difficulty would be significantly exacerbated in the presence of missing data.

Probabilistic models for matrix/tensor factorization have attracted much interest in collaborative filtering and matrix/tensor completion. Probabilistic matrix factorization was proposed in , and its fully Bayesian treatment using Markow chain Monte Carlo (MCMC) inference was shown in and using variational Bayesian inference in . Further extensions of nonparametric and robust variants were presented in . The probabilistic frameworks of tensor factorization were presented in . Other variants include extensions of the exponential family model and the nonparametric Bayesian model . However, the tensor rank or model complexity are often given by a tuning parameter selected by either maximum likelihood or cross-validations, which are computationally expensive and inaccurate. Another important issue is that the inference of factor matrices is performed by either point estimation, which is prone to overfitting, or MCMC inference, which tends to converge very slowly.

To address these issues, we propose a fully Bayesian probabilistic tensor factorization model according to the CP factorization framework. Our objective is to infer the underlying multilinear factors from a noisy incomplete tensor and the predictive distribution of missing entries, while the rank of the true latent tensor can be determined automatically and implicitly. To achieve this, we specify a sparsity-inducing hierarchical prior over multiple factor matrices with individual hyperparameters associated to each latent dimension, such that the number of components in factor matrices is constrained to be minimum. All the model parameters, including noise precision, are considered to be latent variables over which the corresponding priors are placed. Due to complex interactions among multiple factors and fully Bayesian treatment, learning the model is analytically intractable. Thus, we resort to the variational Bayesian inference and derive a deterministic solution to approximate the posteriors of all the model parameters and hyperparameters. Our method is characterized as a tuning parameter-free approach that can effectively avoid parameter selections. The extensive experiments and comparisons on synthetic data illustrate the advantages of our approach in terms of rank determination, predictive capability, and robustness to overfitting. Moreover, several real-word applications, including image completion, restoration, and synthesis, demonstrate that our method outperforms state-of-the-art approaches, including both tensor factorization and tensor completion, in terms of the predictive performance.

The rest of this paper is organized as follows. In Section 2, preliminary multilinear operations and notations are presented. In Section 3, we introduce the probabilistic CP model specification and the model learning via Bayesian inference. A variant of our method using mixture priors is proposed in Section 4. In Section 5, we present the comprehensive experimental results for both synthetic data and real-world applications, followed by our conclusion in Section 6.

Preliminaries and Notations

The inner product of two tensors is defined by $\langle\boldsymbol{\mathcal{A}},\boldsymbol{\mathcal{B}}\rangle=\sum_{i_{1},i_{2},...,i_{N}}\mathcal{A}_{i_{1}i_{2}...i_{N}}\mathcal{B}_{i_{1}i_{2}...i_{N}}$, and the squared Frobenius norm by $\|\boldsymbol{\mathcal{A}}\|_{F}^{2}=\langle\boldsymbol{\mathcal{A}},\boldsymbol{\mathcal{A}}\rangle$. As an extension to $N$ variables, the generalized inner product of a set of vectors, matrices, or tensors is defined as a sum of element-wise products. For example, given $\{\mathbf{A}^{(n)}|n=1,\ldots,N\}$, we define

while the Khatri-Rao product of a set of matrices, except the $n$th matrix, denoted by $\mathbf{A}^{(\backslash n)}$, is

Bayesian Tensor Factorization

Let $\boldsymbol{\mathcal{Y}}$ be an incomplete $N$th-order tensor of size $I_{1}\times I_{2}\times\cdots\times I_{N}$ with missing entries. The element $\mathcal{Y}_{i_{1}i_{2}\ldots i_{N}}$ is observed if $({i_{1},i_{2},\cdots,i_{N}})\in\Omega$, where $\Omega$ denotes a set of indices. For simplicity, we also define a binary tensor $\boldsymbol{\mathcal{O}}$ of the same size as $\boldsymbol{\mathcal{Y}}$ as an indicator of observed entries. We assume $\boldsymbol{\mathcal{Y}}$ is a noisy observation of true latent tensor $\boldsymbol{\mathcal{X}}$, that is, $\boldsymbol{\mathcal{Y}}=\boldsymbol{\mathcal{X}}+\boldsymbol{\mathcal{\varepsilon}}$, where the noise term is assumed to be an i.i.d. Gaussian distribution, i.e., $\boldsymbol{\mathcal{\varepsilon}}\sim\prod_{i_{1},\ldots,i_{N}}\mathcal{N}(0,\tau^{-1})$, and the latent tensor $\boldsymbol{\mathcal{X}}$ can be exactly represented by a CP model, given by

The CP generative model, together with noise assumption, directly give rise to the observation model, which is factorized over observed tensor elements

where the parameter $\tau$ denotes the noise precision, and $\left\langle\mathbf{a}^{(1)}_{i_{1}},\mathbf{a}^{(2)}_{i_{2}},\cdots,\mathbf{a}^{(N)}_{i_{N}}\right\rangle=\sum_{r}\prod_{n}{a}^{(n)}_{i_{n}r}$ denotes a generalized inner-product of $N$ vectors. The likelihood model in (6) indicates that $\mathcal{Y}_{i_{1}\cdots i_{N}}$ is generated by multiple $R$-dimensional latent vectors \big{\{}\mathbf{a}^{(n)}_{i_{n}}\big{|}n=1,\ldots,N\big{\}}, where each latent vector $\mathbf{a}_{i_{n}}^{(n)}$ contributes to a set of observations, i.e., a subtensor whose mode-$n$ index is $i_{n}$. The essential difference between matrix and tensor factorization is that the inner product of $N\geq 3$ vectors allows us to model the multilinear interaction structure, which however leads to many more difficulties in model learning.

In general, the effective dimensionality of the latent space, i.e., ${Rank_{CP}(\boldsymbol{\mathcal{X}})}=R$, is a tuning parameter whose selection is quite challenging and computational costly. Therefore, we seek an elegant automatic model selection, which can not only infer the rank of the latent tensor $\boldsymbol{\mathcal{X}}$, but also effectively avoid overfitting. To achieve this, a set of continuous hyperparameters are employed to control the variance related to each dimensionality of the latent space, respectively. Since the minimum $R$ is desired in the sense of low rank approximation, a sparsity-inducing prior is specified over these hyperparameters, resulting in it being possible to achieve automatic rank determination as a part of the Baybesian inference process. This technique is related to automatic relevance determination (ARD) or sparse Bayesian learning . However, unlike the traditional methods that place the ARD prior over either latent variables or weight parameters, such as Bayesian principle component analysis , our method considers all model parameters as latent variables over which a sparsity-inducing prior is placed with shared hyperparameters.

More specifically, we place a prior distribution over the latent factors, governed by hyperparameters $\boldsymbol{\lambda}=[\lambda_{1},\ldots,\lambda_{R}]$ where each $\lambda_{r}$ controls $r$th component in $\mathbf{A}^{(n)}$, which is

where $\boldsymbol{\Lambda}=\text{diag}(\boldsymbol{\lambda})$ denotes the inverse covariance matrix, also known as the precision matrix, and is shared by latent factor matrices in all modes. We can further define a hyperprior over $\boldsymbol{\lambda}$, which is factorized over latent dimensions

where $\text{Ga}(x|a,b)=\frac{b^{a}x^{a-1}e^{-bx}}{\Gamma(a)}$ denotes a Gamma distribution and $\Gamma(a)$ is the Gamma function.

Since the sparsity is enforced in the latent dimensions, the initialization point of the dimensionality of latent space (i.e., R) is usually set to its maximum possible value, while the effective dimensionality can be inferred automatically under a Bayesian inference framework. It should be noted that since the priors are shared across $N$ latent matrices, our framework can learn the same sparsity pattern for them, yielding the minimum number of rank-one terms. Therefore, our model can effectively infer the rank of tensor while performing the tensor factorization, which can be treated as a Bayesian low-rank tensor factorization.

To complete the model with a fully Bayesian treatment, we also place a hyperprior over the noise precision $\tau$, that is,

For simplicity of notation, all unknowns including latent variables and hyperparameters are collected and denoted together by $\Theta=\{\mathbf{A}^{(1)},\ldots,\mathbf{A}^{(N)},\boldsymbol{\lambda},\tau\}$. The probabilistic graph model is illustrated in Fig. 1, from which we can easily write the joint distribution of the model as

By combining the likelihood in (6), the priors of model parameters in (7), and the hyperpriors in (8) and (9), the logarithm of the joint distribution is given by (see Sec. 1 of Appendix for details)

where $M=\sum_{i_{1},\ldots,i_{N}}\mathcal{O}_{i_{1}\ldots i_{N}}$ denotes the total number of observations. Without loss of generality, we can perform maximum a posteriori (MAP) estimation of $\Theta$ by maximizing (10), which is, to some extent, equivalent to optimizing a squared error function with regularizations imposed on the factor matrices and additional constraints imposed on the regularization parameters.

However, our objective is to develop a method that, in contrast to the point estimation, computes the full posterior distribution of all variables in $\Theta$ given the observed data, that is,

Based on the posterior distribution of $\Theta$, the predictive distribution over missing entries, denoted by $\boldsymbol{\mathcal{Y}}_{\backslash\Omega}$, can be inferred by

2 Model Learning via Bayesian Inference

An exact Bayesian inference in (11) and (12) would integrate over all latent variables as well as hyperparameters, which is obviously analytically intractable. In this section, we describe the development of a deterministic approximate inference under variational Bayesian (VB) framework to learn the probabilistic CP factorization model.

We therefore seek a distribution $q(\Theta)$ to approximate the true posterior distribution $p(\Theta|\boldsymbol{\mathcal{Y}}_{\Omega})$ by minimizing the KL divergence, that is,

where $\ln p(\boldsymbol{\mathcal{Y}}_{\Omega})$ represents the model evidence, and its lower bound is defined by $\mathcal{L}(q)=\int q(\Theta)\ln\left\{\frac{p(\boldsymbol{\mathcal{Y}}_{\Omega},\Theta)}{q(\Theta)}\right\}d\Theta$. Since the model evidence is a constant, the maximum of the lower bound occurs when the KL divergence vanishes, which implies that $q(\Theta)=p(\Theta|\boldsymbol{\mathcal{Y}}_{\Omega})$.

For the initial derivation, it will be assumed that the variational distribution is factorized w.r.t. each variable $\Theta_{j}$ and therefore can be written as

It should be noted that this is the only assumption about the distribution, while the particular functional forms of the individual factors $q_{j}(\Theta_{j})$ can be explicitly derived in turn. The optimised form of the $j$th factor based on the maximization of $\mathcal{L}(q)$ is given by

As can be seen from the graphical model shown in Fig. 1, the inference of mode-$n$ factor matrix $\mathbf{A}^{(n)}$ can be performed by receiving the messages from observed data and its co-parents, including other factors $\mathbf{A}^{(k)},k\neq n$ and the hyperparameter $\tau$, which are expressed by the likelihood term (6), and incorporating the messages from its parents, which are expressed by the prior term (7). By applying (15), it has been shown that their posteriors can be factorized as independent distributions of their rows, which are also Gaussian (see Sec. 2 of Appendix for details), given by

where the posterior parameters can be updated by

For simplicity, we attempt to compute (19) by multilinear operations. Let $\forall n$, $\mathbf{B}^{(n)}$ of size $I_{n}\times R^{2}$ denote an expectation of a quadratic form related to $\mathbf{A}^{(n)}$ by defining the $i_{n}$th-row vector as

where $\mathbf{1}_{\prod_{n}{I_{n}}}$ denotes a vector of length $\prod_{n}I_{n}$ and all elements are equal to one.

where $\boldsymbol{\mathcal{O}}_{\cdots i_{n}\cdots}$ denotes a subtensor by fixing model-$n$ index to $i_{n}$. It should be noted that the Khatri-Rao product is computed by all mode factors except the $n$th mode, while the sum is performed according to the indices of observations, implying that only factors that interact with $\mathbf{a}^{(n)}_{i_{n}}$ are taken into account. Another complicated part in (17) can also be simplified by multilinear operations, i.e.,

2.2 Posterior distribution of hyperparameters 𝝀𝝀\boldsymbol{\lambda}

It should be noted that, instead of point estimation via optimizations, learning the posterior of $\boldsymbol{\lambda}$ is crucial for automatic rank determination. As seen in Fig. 1, the inference of $\boldsymbol{\lambda}$ can be performed by receiving messages from $N$ factor matrices and incorporating the messages from its hyperprior. By applying (15), we can identify the posteriors of $\lambda_{r},\forall r\in[1,R]$ as an independent Gamma distribution (see Sec. 4 of Appendix for details),

where $c_{M}^{r}$, $d_{M}^{r}$ denote the posterior parameters learned from $M$ observations and can be updated by

The expectation of the inner product of the $r$th component in mode-$n$ matrix w.r.t. $q$ distribution can be evaluated using the posterior parameters in (16), i.e.,

By combining (25) and (26), we can further simplify the computation of $\mathbf{d}_{M}=[d_{M}^{1},\ldots d_{M}^{R}]^{T}$ as

2.3 Posterior distribution of hyperparameter τ𝜏\tau

The inference of the noise precision $\tau$ can be performed by receiving the messages from observed data and its co-parents, including $N$ factor matrices, and incorporating the messages from its hyperprior. By applying (15), the variational posterior is a Gamma distribution (see Sec. 5 of Appendix for details), given by

where the posterior parameters can be updated by

However, the posterior expectation of model error in the above expression cannot be computed straightforwardly, and therefore, we need to introduce the following results.

Assume a set of independent $R$-dimensional random vectors $\{\mathbf{x}^{(n)}|n=1,\ldots,N\}$, then

where the left term denotes the expectation of the squared inner product of $N$ vectors, and the right term denotes the inner product of $N$ matrices, where each matrix of size $R\times R$ denotes an expectation of the outer product of the $n$th vector, respectively.

Given a set of independent random matrices $\{\mathbf{A}^{(n)}|n=1,\ldots,N\}$, we assume that $\forall n,\forall i_{n}$, the row vectors $\{\mathbf{a}_{i_{n}}^{(n)}\}$ are independent, then

From Theorems 3.2 and 3.3, the posterior expectation term in (29) can be evaluated explicitly. Due to the missing entries in $\boldsymbol{\mathcal{Y}}$, the evaluation form is finally written as (see Sec. 8 of Appendix for details)

An intuitive interpretation of (29) is straightforward. $a_{M}$ is related to the number of observations and $b_{M}$ is related to the residual of model fitting measured by the squared Frobenius norm on observed entries.

2.4 Lower bound of model evidence

The inference framework presented in the previous section can essentially maximize the lower bound of model evidence that is defined in (13). Since the lower bound should not decrease at each iteration, it can be used to test for convergence. The lower bound of the log-marginal likelihood is computed by

where the first term denotes the posterior expectation of joint distribution, and the second term denotes the entropy of posterior distributions.

Various terms in the lower bound are evaluated and derived by taking parametric forms of $q$ distribution, giving the following results (see Sec. 9 of Appendix for details)

An intuitive interpretation of (35) is as follows. The first term is related to model residual; the second term is related to the weighted sum of squared $L_{2}$-norm of each component in factor matrices, while the uncertainty information is also considered; the rest terms are related to negative KL divergence between the posterior and prior distributions of hyperparameters.

2.5 Initialization of model parameters

2.6 Interpretaion of automatic rank determination

The entire procedure of model inference is summarized in Algorithm 1. It should be noted that tensor rank is determined automatically and implicitly. More specifically, updating $\boldsymbol{\lambda}$ in each iteration results in a new prior over $\{\mathbf{A}^{(n)}\}$, and then, $\{\mathbf{A}^{(n)}\}$ can be updated using this new prior in the subsequent iteration, which in turn affects $\boldsymbol{\lambda}$. Hence, if the posterior mean of $\lambda_{r}$ becomes very large, the $r$th components in $\{\mathbf{A}^{(n)}\},\forall n\in[1,N]$ are forced to be zero because of their prior information, and the tensor rank can be obtained by simply counting the number of non-zero components in the factor matrices. For implementation of the algorithm, we can keep the size of $\{\mathbf{A}^{(n)}\}$ unchanged during iterations; an alternative method is to eliminate the zero-components of $\{\mathbf{A}^{(n)}\}$ after each iteration.

3 Predictive Distribution

The predictive distributions over missing entries, given observed entries, can be approximated by using variational posterior distribution, that is,

Thus, the predictive variance can be obtained by $\text{Var}(\mathcal{Y}_{i_{1}\ldots i_{N}})=\frac{\nu_{y}}{\nu_{y}-2}\mathcal{S}_{i_{1}\ldots i_{N}}^{-1}$.

4 Computational Complexity

The computation cost of the $N$ factor matrices in (17) is $O(NR^{2}M+R^{3}\sum_{n}I_{n})$, where $N$ is the order of the tensor, $M$ denotes the number of observations, i.e., the input data size. $R$ is the number of latent components in each $\mathbf{A}^{(n)}$, i.e., model complexity or tensor rank, and is generally much smaller than the data size, i.e., $R\ll M$. Hence, it has linear complexity w.r.t. the data size and polynomial complexity w.r.t. the model complexity. It should be noted that, because of the automatic model selection, the excessive latent components are pruned out in the first few iterations such that $R$ reduces rapidly in practice. The computation cost of the hyperparameter $\boldsymbol{\lambda}$ in (25) is $O(R^{2}\sum_{n}I_{n})$, which is dominated by the model complexity, while the computation cost of noise precision $\tau$ in (29) is $O(R^{2}M)$. Therefore, the overall complexity of our algorithm is $O(NR^{2}M+R^{3})$, which scales linearly with the data size but polynomially with the model complexity.

5 Discussion of Advantages

The advantages of our method are discussed as follows:

The automatic determination of CP rank enables us to obtain an optimal low-rank tensor approximation, even from a highly noisy and incomplete tensor.

Our method is characterized as a tuning parameter-free approach and all model parameters can be inferred from the observed data, which avoids the computational expensive parameter selection procedure. In contrast, the existing tensor factorization methods require a predefined rank, while the tensor completion methods based on nuclear norm require several tuning parameters.

The uncertainty information over both latent factors and predictions of missing entries can be inferred by our method, while most existing tensor factorization and completion methods provide only the point estimations.

An efficient and deterministic Bayesian inference is developed for model learning, which empirically shows a fast convergence.

Mixture Factor Priors

The low-rank assumption is powerful in general cases, however if the tensor data does not satisfy an intrinsic low-rank structure and a large amount of entries are missing, it may yield an oversimplified model. In this section, we present a variant of Bayesian CP factorization model which can take into account the local similarity in addition to the low-rank assumption.

We specify a Gaussian mixture prior over factor matrices such that the prior distribution in (7) can be rewritten as $\forall i_{n}\in[1,I_{n}],\forall n\in[1,N]$,

Experimental Results

We conducted extensive experiments using both synthetic data and real-world applications, and compared our fully Bayesian CP factorization (FBCP)Matlab codes are available at http://www.bsp.brain.riken.jp/~qibin/homepage/BayesTensorFactorization.html with seven state-of-the-art methods. Tensor factorization based scheme includes CPWOPT and CPNLS , while the completion based scheme includes HaLRTC and FaLRTC , FCSA , hard-completion (HardC.) , geomCG and STDC . Our objective when using synthetic data was to validate our method from several aspects: i) capability of rank determination; ii) reconstruction performance given a complete tensor; iii) predictive performance over missing entries given an incomplete tensor. Two real-world applications including image inpainting and facial image synthesis were used for demonstration. All experiments were performed by a PC (Intel Xeon(R) 3.3GHz, 64GB memory).

The synthetic tensor data is generated by the following procedure. $N$ factor matrices $\{\mathbf{A}^{(n)}\}_{n=1}^{N}$ are drawn from a standard normal distribution, i.e., $\forall n,\forall i_{n},\mathbf{a}^{(n)}_{i_{n}}\sim\mathcal{N}(\mathbf{0},\mathbf{I}_{R})$. Then, the true tensor is constructed by $\boldsymbol{\mathcal{X}}=[\![\mathbf{A}^{1},\ldots,\mathbf{A}^{(N)}]\!]$, and an observed tensor by $\boldsymbol{\mathcal{Y}}=\boldsymbol{\mathcal{X}}+\boldsymbol{\mathcal{\varepsilon}}$, where $\boldsymbol{\mathcal{\varepsilon}}\sim\prod_{i_{1},\ldots,i_{N}}\mathcal{N}(0,\sigma^{2}_{\boldsymbol{\mathcal{\varepsilon}}})$ denotes i.i.d. additive noise. The missing entries, chosen uniformly, are marked by an indicator tensor $\boldsymbol{\mathcal{O}}$.

To illustrate our model, we provide two demo videos in the supplemental materials. A true latent tensor $\boldsymbol{\mathcal{X}}$ is of size $10\times 10\times 10$ with CP rank $R=5$, the noise parameter was $\sigma^{2}_{\boldsymbol{\mathcal{\varepsilon}}}=0.001$, and $40\%$ of entries were missing. Then, we applied our method with the initial rank being set to 10. As shown in Fig. 2, three factor matrices are inferred in which five components are effectively pruned out, resulting in correct estimation of tensor rank. The lower bound of model evidence increases monotonically, which indicates the effectiveness and convergence of our algorithm. Finally, the posterior of noise precision $\tau\approx 1000$ implies the method’s capability of denoising and the estimation of $\sigma^{2}_{\boldsymbol{\mathcal{\varepsilon}}}\approx 0.001$, $\text{SNR}=10\log\frac{\sigma^{2}_{\boldsymbol{\mathcal{X}}}}{\tau^{-1}}$.

1.2 Automatic determination of tensor rank

To evaluate the automatic determination of tensor rank (i.e., CP rank), extensive simulations were performed under varying experimental conditions related to tensor size, tensor rank, noise level, missing ratio, and the initialization method of factor matrices (e.g., SVD or random sample). Each result is evaluated by 50 repetitions corresponding to 50 different tensors generated under the same criterion. There are four groups of experiments. (A) Given complete tensors of size $20\times 20\times 20$ with $R=5$, the evaluations were performed under varying noise levels and by two different initializations (see Fig. 3). (B) Given incomplete tensors of size $20\times 20\times 20$ with $R=5$ and SNR=20 dB, the evaluations were performed under five different missing ratios, and by different initializations (see Fig. 3). (C) Given incomplete tensors with $R=5$ and SNR=0 dB, the evaluations were performed under varying missing ratios and two different tensor sizes (see Fig. 3). (D) Given incomplete tensors of size $20\times 20\times 20$ with SNR=20 dB, the evaluations were performed under varying missing ratios and two different true ranks (see Fig. 3).

From Fig. 3, we observe that SVD initialization is slightly better than random initialization in terms of the determination of tensor rank. If the tensor is complete, our model can detect the true tensor rank with $100\%$ accuracy when SNR$\geq$10 dB. Although the accuracy decreased to 70% under a high noise level of 0 dB, the error deviation is only $\pm 1$. On the other hand, if the tensor is incomplete and almost free of noise, the detection rate is 100%, when missing ratio is 0.7, and is 44% with an error deviation of only $\pm 1$, even under a high missing ratio of 0.9. As both missing data and high noise level are presented, our model can achieve 90% accuracy under the condition of SNR=0 dB and 0.5 missing ratio. It should be noted that, when the data size is larger, such as $50\times 50\times 50$, our model can achieve 90% accuracy, even when SNR=0 dB and the missing ratio is 0.9. If the true rank is larger, such as $R=15$, the model can correctly recover the rank from a complete tensor, but fails to do so when the missing ratio is larger than 0.5.

We can conclude from these results that the determination of the tensor rank depends primarily on the number of observed entries and the true tensor rank. In general, more observations are necessary if the tensor rank is larger; however, when high-level noise occurs, the excessive number of observations may not be helpful for rank determination.

1.3 Predictive performance

In this experiment, we considered incomplete tensors of size $20\times 20\times 20$ generated by the true rank $R=5$ and SNR=30 dB under varying missing ratios. The initial rank was set to 10. The relative standard error $RSE=\frac{\|\hat{\boldsymbol{\mathcal{X}}}-\boldsymbol{\mathcal{X}}\|_{F}}{\|\boldsymbol{\mathcal{X}}\|_{F}}$, where $\hat{\boldsymbol{\mathcal{X}}}$ denotes the estimation of the true tensor $\boldsymbol{\mathcal{X}}$, was used to evaluate the performance. To ensure statistically consistent results, the performance is evaluated by 50 repetitions for each condition. As shown in Fig. 4, our method significantly outperforms other algorithms under all missing ratios. Factorization-based methods, including CPWOPT, and CPNLS show a better performance than completion-based methods when the missing ratio is relatively small, while they perform worse than completion methods when the missing ratio is large, e.g., 0.9. FaLRTC, FCSA, and HardC. achieve similar performances, because they are all based on nuclear norm optimization. geomCG achieves a performance comparable with that of CWOPT and CPNLS when data is complete, while it fails as the missing ratio becomes high. This is because geomCG requires a large number of observations and precisely defined rank. It should be noted that STDC outperforms all algorithms except FBCP as the missing ratio becomes extremely high. These results demonstrate that FBCP, as a tensor factorization method, can be also effective for tensor completion, even when an extremely sparse tensor is presented.

We also conducted two additional experiments. One is the reconstruction from a complete tensor, the other is the tensor completion when the noise level is high, i.e., SNR =0 dB. The results of these two experiments are presented in Appendix (see Sec. 11, 12).

2 Image Inpainting

In this section, the applications of image inpainting based on several benchmark images, shown in Fig. 5, are used to evaluate and compare the performance of different methods. The colorful image can be represented by a third-order tensor of size $200\times 200\times 3$. We conducted various experiments under four groups of conditions. (A) Structural image with uniformly random missing pixels. A building facade image with 95% missing pixels under two noise conditions, i.e., noise free and SNR=5dB, were considered as observations. (B) Natural image with uniformly random missing pixels. The Lena image of size $300\times 300$ with 90% missing pixels under two noise conditions, i.e., noise free and SNR=10dB, were considered. (C) Non-random missing pixels. We conducted two experiments for image restoration from a corrupted image: 1) The Lenna image corrupted by superimposed text was used as an observed imageA demo video is available in the supplemental materials.. In practice, the location of text pixels are difficult to detect exactly; we can simply indicate missing entries by a value larger than 200 to ensure that the text pixels are completely missing. 2) The scrabbled Lenna image was used as an observed image and pixels with values larger than 200 can be marked as missing. (D) Object removal. Given an image and a mask covering the object area, our goal was to complete the image without that object. The algorithm settings of compared methods are described as follows. For factorization-based methods, the initial rank was set to 50 in cases of (A) and (B) due to the high missing ratios, and 100 in cases of (C) and (D). For completion-based methods, the tuning parameters were chosen by multiple runs and performance evaluated on the ground-truth of missing pixels.

The visual effects of image inpainting are shown in Fig. 6, and the predictive performances are shown in Table I where case (D) is not available due to the lack of ground-truth. In case (A), we observe that FBCP outperforms other methods for a structural image under an extremely high missing ratio and the superiority is more significant when an additive noise is involved. In case (B), observe that STDC obtains the best performance followed by FBCP that is better than other methods. However, STDC is severely degraded when noise is involved, and obtains the same performance as FBCP, while its visual quality is still much better than others. These indicate that the additional smooth constraints in STDC are suitable for natural image. In case (C), FBCP is superior to all other methods, followed by STDC. The completion-based methods obtain relatively smoother effects than factorization-based methods, but the global color of the image is not recovered well, resulting in a poor predictive performance. In case (D), FBCP obtains the most clean image by removing the object completely while the ghost effects appear in all other methods. HaLRTC, FaLRTC and FCSA outperform CPWOPT, CPNLS and STDC.

From these results we can conclude that the completion-based methods generally outperforms factorization-based methods for image completion. However, FBCP significantly improves ability of factorization-based scheme by automatic model selection and robustness to overfitting, resulting in potential applications for various image inpainting problems. The necessary number of observed entries mainly depends on the rank of true image. For instance, a structural image with an intrinsic low-rank need very fewer observations than a natural image. However only 10% observed pixels from lena image are not sufficient to recover the whole image, which is caused by the low-rank assumption. This property is common for all these algorithms except STDC, because STDC employs an auxiliary information as additional constraints. The advantage of STDC has been shown for lena image, while its disadvantages are that the auxiliary information must be well designed for each specific application, which make it difficult to be applied to other types of data. In addition, STDC degrades in presence of the non-random missing pixels or noise. Moreover, the performance of HaLRTC and STDC are sensitive to tuning parameters that must be carefully selected for each specific condition. Therefore, a crucial drawback of completion-based scheme lies in the tuning parameters whose selection is quite challenging when the ground-truth of missing data is unknown.

Next, we perform image completion extensively on eight images in Fig. 5 with randomly missing pixels. Since most of these images are natural images on which the low-rank approximation cannot recover the missing pixels well, we apply the fully Bayesian CP with mixture priors (FBCP-MP) for comparison with FBCP and other related methods. For FBCP and FBCP-MP, the same initialization of $R=100$ was applied, while CPWOPT was performed by using the optimal ranks obtained from FBCP and FBCP-MP to show the best performance. The parameter selection for other methods was same with previous experiments. The size of all images is $256\times 256\times 3$. Table II shows quantitative results in terms of recovery performance and runtime. Observe that FBCP-MP improves the performance of FBCP significantly and achieves the best recovery performance, especially in the case of high missing rate. The time costs of FBCP and FBCP-MP are comparable with completion-based methods and significantly lower than other tensor factorization method. STDC obtains the comparable performance with FBCP-MP, however the parameters must be manually tuned for the specific condition. More detailed results on each image are shown visually and quantitatively in the supplemental materials. These results demonstrate the effectiveness of mixture priors and advantages when the local similarity is taken into account in addition to the low-rank assumption.

3 Facial Image Synthesis

For recognition of face images captured from surveillance videos, the ideal solution is to create a robust classifier that is invariant to some factors, such as pose and illumination. Hence, there arises the question whether we can generate novel facial images under multiple conditions given images under other conditions. Tensors are highly suitable for modeling a multifactor image ensemble, and therefore, we introduce a novel application of facial image synthesis that utilizes tensor factorization approaches.

We used the dataset of 3D Basel Face Model , which contains an ensemble of facial images of 10 people, each rendered in 9 different poses under 3 different illuminations. All 270 facial images were decimated and cropped to $68\times 68$ pixels, and were then represented by a fourth-order tensor of size $4624\times 10\times 9\times 3$. As shown in Fig. 7, some images were fully missing. Since some methods are either computationally intractable or not applicable to $N\geq 4$ order tensor, five algorithms were finally applied on this dataset under different missing ratios. The initial rank was set to 100 in factorization based methods, while the parameters of completion based methods were well tuned based on the ground-truth of missing images.

As shown in Fig. 8, the visual effects of image synthesis produced by FBCP are significantly superior to those produced by other methods. Although both CPWOPT and HaLRTC produce images that are too smooth and blurred, HaLRTC obtains much better visual quality than CPWOPT. The detailed performances are compared in Table III, where RSE w.r.t. observed entries reflects the performance of model fitting, and RSE w.r.t. missing entries particularly reflects the predictive ability. Note that RSE $=$N/A implies that HaLRTC and FaLRTC donot model the observed entries. The inferred rank by FBCP is within the range of that are close to the initialization. Observe that completion based methods including HaLRTC, FaLRTC and HardC. achieve better performance than CPWOPT. However, FBCP demonstrates the possibility that factorization-based scheme can significantly outperform completion-based methods, especially in terms of performance on missing images.

Conclusion

In this paper, we proposed a fully Bayesian CP factorization which can naturally handle incomplete and noisy tensor data. By employing hierarchical priors over all unknown parameters, we derived a deterministic model inference under a fully Bayesian treatment. The most significant advantage is automatic determination of CP rank. Moreover, as a tuning parameter-free approach, our method avoids the parameter selection problem and can also effectively prevent overfitting. In addition, we proposed a variant of our method by using mixture priors, which shows advantages on natural images with a large amount of missing pixels. Empirical results validate the effectiveness in terms of discovering the ground-truth of tensor rank and imputing missing values for an extremely sparse tensor. Several real-world applications, such as image completion and image synthesis, demonstrated the superiority of our method over state-of-the-art techniques. Due to several interesting properties, our method would be attractive for many potential applications.

References