In this paper we consider the tensor completion problem. We suppose that values of tensor are generated by some smooth function, i.e.
where is a point on some multi-dimensional grid and is some unknown smooth function. However, the tensor values are known only at some small subset of the grid. The task is to complete the tensor, i.e., to reconstruct the tensor values at all points on the grid taking into account the properties of the data generating process .
This problem statement differs from the traditional problem statement, which does not use any assumptions about the function . Knowing some properties of the data generating function provides insights about how the tensor values relate to each other, and this, in turn, allows us to improve the results. In this work we assume that function is smooth.
There are a lot of practical applications that suit the statement. For example, modeling of physical processes, solutions of differential equations, modeling probability distributions.
In this paper we propose to model the smoothness of the data generating process by using Gaussian Process Regression (GPR). In GPR the assumptions about the function that we approximate are controlled via the kernel function. The GPR model is then used to construct the initial solution to the tensor completion problem
In principle, such initialization can improve any other tensor completion technique. It means that using the proposed initialization state-of-the-art results can be obtained employing some simple optimization procedure like Stochastic Gradient Descent.
When the tensor order is high the problem should be solved in some low-rank format because the number of elements of the tensor grows exponentially. The proposed approach is based on the tensor train (TT) format for its computational efficiency and ability to handle large dimensions .
The contributions of this paper are as follows
We introduce new initialization algorithm which takes into account the tensor generating process. The proposed algorithm is described in Section 3.
The proposed initialization technique automatically selects the rank of the tensor, the details are given in Section 3.2.
We conducted empirical evaluations of the proposed approach and compared it with tensor completion techniques without our initialization. The results are given in Section 4 and show the superiority of the proposed algorithm.
2 Tensor completion
The formal problem statement is as follows. Suppose that is a -way tensor, (by tensor here we mean a multi-dimensional array). Tensor values are known only at some subset of indices . By we denote the projection onto the set , i.e.
We formulate the tensor completion as an optimization problem
where is a tensor train rank of , which is a generalization of the matrix rank, and is the Frobenius norm. A tensor is said to be in tensor train format if its elements are represented as
where is a three-way tensor core with size ,
. Vectoris called TT-rank.
Tensor train format assumes that the full tensor can be approximated by a set of -way core tensors, the total number of elements in core tensors is , where , , which is much smaller than .
In problem (1) we optimize the objective function straightforwardly with respect to tensor cores while having their sizes fixed. Problem (1) is non-convex, so optimization methods can converge to a local minimum. To get an efficient solution we impose two requirements:
Initial tensor in tensor train format should be as close to the optimum as possible.
Availability of an efficient optimization procedure that will be launched from the obtained initial tensor.
These steps are independent, and one can apply any desired algorithm in each of them.
In this work we develop the initialization algorithm that allows obtaining accurate initial tensor for the case when the tensor of interest is generated by some smooth function. The experimental section demonstrates that our initialization can improve the results of any optimization procedure, ensuring our approach to be universal.
We consider tensors that are generated by some function, i.e. tensor values are computed as follows
where is some unknown smooth function and , , are tensor sizes. The set of points is a full factorial Design of Experiments, i.e. a multi-dimensional grid, and we also assume that the grid is uniform.
In this setting the tensor completion can be considered as a regression problem and can be solved by any regression technique that guarantees the smoothness of the solution. However, in the tensor completion problem we are interested in a tensor of values of at a predefined finite grid of points. The tensor should be in a low-rank format to be able to perform various operations with tensor efficiently (e.g. calculation of the norm of the tensor, dot product and other).
These observations give us the solution — build regression model using the observed values of the tensor, then use the obtained approximation as a black-box for the TT-cross approximation algorithm . The last step results in a tensor in TT format, which is a low-rank format and allows efficient computations. The next step (which is optional) is to improve the obtained solution by using it as initialization for any other tensor completion technique.
Let us write down the set of observed tensor values into a vector and the corresponding indices into a matrix (each row is a vector of indices ). Then the approach for tensor completion (in TT format) can be written as follows
Apply some tensor completion technique using as an initial value.
At step (a) the choice of the regression technique affects the result of the initialization, although it can be arbitrary. It is required to choose the regression algorithm such that it will capture the peculiarities of the tensor we would like to restore. In this work we suppose that the tensor generating function is smooth (which is a common situation when modeling physical processes). Therefore, we choose a regression technique that is good at approximating smooth functions. A reasonable choice, in this case, is to use Gaussian Process (GP) Regression . GP models is a favorite tool in many engineering applications as they proved to be efficient, especially for problems where it is required to model some smooth function . The points are not given, all we know is that at the point with multi-index on the grid the function value is equal to . To make the problem statement reasonable we assume that the indices are connected with the points as follows: , where . So, as an input for the approximation we set and such that .
At step (b) we use TT-cross because it allows to efficiently approximate black-box function by a low-rank tensor in TT format. Moreover, this approach can automatically select TT-rank making it more desirable. More details on the technique are given in Section 3.2.
The described approach has the following benefits
Initial tensor which is close to the optimal value in terms of the reconstruction error at observed values. It will push the optimization to faster convergence.
Better generalization ability — there are many degrees of freedom: a lot of different tensor train factors can give low reconstruction error at observed positions but can give a large error at other locations. Accurate approximation model will push the initial tensor to be closer to the original tensor in both the observed positions and unobserved ones.
TT-cross technique chooses rank automatically, so there is no need to tune the rank of the tensor manually.
3.1 Gaussian Process Regression
One of the most efficient tools for approximating smooth functions is the Gaussian Process (GP) Regression . GP regression is a Bayesian approach where a prior distribution over continuous functions is assumed to be a Gaussian Process, i.e.
where is a vector of outputs, is a matrix of inputs, ,
is a noise variance,is a mean vector modeled by some function , is a covariance matrix for some a priori selected covariance function and
is an identity matrix. An example of such function is a squared exponential kernel
are parameters of the kernel (hyperparameters of the GP model). The hyperparameters should be chosen according to the given data set.
Without loss of generality we make the standard assumption of zero-mean data. Now, for a new unseen data point we have
where and .
Let us denote the vector of hyperparameters and by . To choose the hyperparameters of our model we consider the log-likelihood
and optimize it over the hyperparameters . The runtime complexity of learning GP regression is as we need to calculate the inverse of , its determinant and derivatives of the log-likelihood. If the sample size ( in our case) is large, the computational complexity becomes an issue. There are several ways to overcome it. If the data set has a factorial structure (multidimensional grid in a simple case) we can use the algorithm from . If the structure is factorial with a small number of missing values the method from  should be applied. For general unstructured cases, the approximate GP model can be built using, for example, the model described in  or use a subsample as a training set.
After tuning of the hyperparameters, we can use the mean of the posterior distribution (2) as a prediction of the model.
Note, that the input points in our case is a set of indices of the observed values . For the GP model we scale each index to interval.
3.2 Tensor-Train cross-approximation
To approximate tensor generated by we use Tensor-Train cross-approximation. First, let us consider the matrix case. Suppose that we are given a rank- matrix of size . A cross-approximation for the matrix is represented as
where are some columns and rows of the matrix and is the submatrix on the intersection of these rows and columns. To construct accurate approximation it is required to find submatrix of large volume. It can be done in operations .
Now for tensor the procedure is the following. At the first step let us consider unfolding of size and rank . Using row-column alternating algorithm from  we can find linearly independent columns of matrix , these columns form matrix . After that applying maxvol procedure  to the matrix we can find set of row indices , matrix and matrix that will give the cross-approximation of unfolding :
where is of size . Next, let us form tensor from rows of :
and reshape it into a tensor of size . Next step is to apply the same procedure to the unfolding of the tensor and obtain the matrices , and
of size .
Repeating the described procedure times we will end up with matrices of sizes . Then each matrix can be reshaped to the -way tensor of size , and can be used as core tensors for TT format. It turns out that such representation is a TT decomposition of the initial tensor .
The exact ranks
are not known to us in general. They can only be estimated from the above (e.g., by the maximum rank of the corresponding unfolding). If the rank is overestimated then the calculation of matricesbecomes unstable operation (because we obtain almost rank-deficient unfolding matrices). However, in  the authors suggest some simple modification that overcomes this issue. Therefore, we need to estimate the ranks from the above, but the estimate should not be much larger than the real rank. So, the approach is to start from some small rank, construct the tensor in TT format and then apply recompression (see ). If there is a rank that is not reduced, then we underestimated that rank and should increase it and repeat the procedure.
4 Experimental results
In this section we present the results of the application of our approach to engineering problems.
The experimental setup is the following. We try the following optimization algorithms
SGD – stochastic gradient descent ,
Ropt – Riemannian optimization ,
TTWopt – weighted tensor train optimization ,
ALS – alternating least squares .
We run each algorithm with random initialization and with the proposed GP-based initialization and then compare results.
The quality of the methods is measured using mean squared error (MSE)
where is some set of indices independent from the given observed set of indices , is a size of the set and is an obtained approximation of the actual tensor .
The listed above algorithms were compared on two problems: CMOS oscillator model and Cookie problem (see Section 4.2 and Section 4.1 correspondingly). In CMOS oscillator problem we run each optimization
times with different random training sets and then calculate the average reconstruction error as well as standard deviation. However, for Cookie problem we performedruns, and the training set was the same during all runs because it is more computationally intensive and generating several training sets takes more time.
Note, that when we use GP based initialization, the TT-rank of the tensor is selected automatically by the TT-cross algorithm and max value of can be larger than . The optimization algorithms with random initialization do not have a procedure for automatic rank selection, so we ran them with different ranks (from to ) and then chose the best one.
TTWopt implementation111https://github.com/yuanlonghao/T3C_tensor_completion does not support high-dimensional problems. For higher dimensional problems the authors of TTWopt propose to use SGD. The authors of TTWopt also propose truncated SVD based initialization. The idea is to fill missing values using the mean value of the observed part of the tensor and then to apply truncated SVD to obtain TT cores. However, such approach is only applicable to low-dimensional tensors as it requires to calculate full matrices of large size.
For Ropt and ALS we used publically available MATLAB codes 222https://anchp.epfl.ch/index-html/software/ttemps/.
4.1 Cookie problem
is a disk of radius and is a number of disks which form grid. This is a heat equation where heat conductivity depends on (see illustration in Fig. 1) and is an -dimensional.
We are interested in average temperature over : . If takes possible values then there are possible values of .
In this work we used the following setup for the Cookie problem: each parameter lie in the interval , number of levels for each p is , number of cookies is and , size of the observed set is , for the test set we used independently generated points.
The results of tensor completion are presented in Table 1. One can see that GP based initialization gives lower reconstruction errors both on the training set and test set except for ALS technique. ALS method with the proposed initialization overfits: the error on the training set is close to , whereas the test error is much more significant. The error on the training set is about , which means that the training set was approximated with machine precision. It is not surprising if we recall that there are only observed values, while the number of free parameters that are used to construct TT is much higher.
4.2 CMOS ring oscillator
Let us consider the CMOS ring oscillator . It is an electric circuit which consists of stages of CMOS inverters. We are interested in the oscillation frequency of the oscillator. The characteristics of the electric circuit are described by parameters. Each parameter can take one of values, so the total size of the tensor is . The number of observed values that were used during the experiments is . For the test set we used independently generated points.
The results of the experiments are given in Table 2. The table demonstrates that utilizing GP based initialization improves the results for all algorithms except ALS. ALS, in this case, overfits again: training error is extremely small, whereas the test error is much larger, though it is rather small compared to other techniques and ALS with random initialization.
|Random init||GP init|
|error||average N iters||error||average N iters|
|Random init||GP init|
|error||N iters||error||N iters|
All in all, the obtained results prove that GP based initialization allows improving the tensor completion results in general. At least it provides better training error. As for the error on the test set one should be more careful as the number of degrees of freedom is large and there exist many solutions that give a small error for the observed values but large errors for other values.
5 Related works
One set of approaches to tensor completion is based on nuclear norm minimization. The nuclear norm of a matrix is defined as a sum of all singular values of the matrix. This objective function is a convex envelope of the rank function. For a tensor the nuclear norm is defined as a sum of singular values of matricizations of the tensor.
There are efficient off-the-shelf techniques for such types of problems that apply interior-point methods. However, they are second-order methods and scale poorly with the dimensionality of the problem. Special optimization technique was derived for nuclear norm minimization [8, 13, 18].
More often such techniques are applied to matrices or low-dimensional tensors as their straightforward formulation allows finding the full tensor. It becomes infeasible when we come to high-dimensional problems.
The second type of approaches is based on low-rank tensor decomposition [1, 7, 12, 19, 23]. There are several tensor decompositions, and all these papers derive some optimization procedure for one of them, namely, CP decomposition, Tucker decomposition or TT/MPS decomposition. The simplest technique is the alternating least squares . It just finds the solution iteratively at each iteration minimizing the objective function w.r.t. one core while other cores are fixed.
Another approach is based on Riemannian optimization that tries to find the optimal solution on the manifold of low-rank tensors of the given structure . The same can be done by using Stochastic Gradient Descent . Riemannian optimization, TTWopt, ALS and its modifications (e.g., ADF, alternating directions fitting ) try to find the TT representation of the actual tensor iteratively. At each iteration it optimizes TT cores such that the resulting tensor approximates well the tensor which coincides with the real tensor at observed indices and with the result of the previous iteration at other indices. All these approaches need to specify rank manually. The authors of  apply the Bayesian framework for CP decomposition which allows them to select the rank of the decomposition automatically.
In some papers the objective is modified by introducing special regularizers to suit the problem better . For example, in [7, 26] to obtain better results for visual data a special prior regularizer was utilized.
Our proposed algorithm is an initialization technique for the tensor completion problems in TT format and can be used with most of the algorithms solving such problems. If the assumptions from Section 3 (the tensor values are values of some rather smooth function of tensor indices) are satisfied the initial value will be close to the optimal providing better results. The question of a good initialization is rarely taken into account. In paper 
a special initialization if proposed for visual data. The idea is to use some crude technique (like bilinear interpolation) to fill missing values and after that apply SVD-based tensor train decomposition. The drawback of the approach is that it can be applied only in case of small-dimensional tensors as we need to fill all missing values. In they propose special initialization for the Alternating Direction Fitting (ADF) method. This is a general technique for the tensor completion and it does not take into account the assumptions on the data generating function.
We proposed a new initialization algorithm for high-dimensional tensor completion in TT format. The approach is designed mostly for the cases when some smooth function generates the tensor values. It can be combined with any optimization procedure that is used for tensor completion. Additionally, the TT-rank of the initial tensor is adjusted automatically by the TT-cross method and defines the resulting rank of the tensor. So, the approach provides an automatic rank selection. Our experimental study confirms that the proposed initialization delivers lower reconstruction errors for many of the optimization procedures.
We would like to thank Zheng Zhang, who kindly provided us CMOS ring oscillator data set.
-  E. Acar, D. M. Dunlavy, T. G. Kolda, and M. Mørup, Scalable tensor factorizations for incomplete data, Chemometrics and Intelligent Laboratory Systems, 106 (2011), pp. 41–56.
-  J. Ballani and L. Grasedyck, Hierarchical tensor approximation of output quantities of parameter-dependent pdes, SIAM/ASA Journal on Uncertainty Quantification, 3 (2015), pp. 852–872.
-  M. Belyaev, E. Burnaev, E. Kapushev, M. Panov, P. Prikhodko, D. Vetrov, and D. Yarotsky, Gtapprox: Surrogate modeling for industrial design, Advances in Engineering Software, 102 (2016), pp. 29–39.
M. Belyaev, E. Burnaev, and Y. Kapushev, Gaussian process regression
for structured data sets
, in International Symposium on Statistical Learning and Data Sciences, Springer, 2015, pp. 106–115.
-  M. Belyaev, E. Burnaev, and Y. Kapushev, Computationally efficient algorithm for gaussian process regression in case of structured samples, Computational Mathematics and Mathematical Physics, 56 (2016), pp. 499–513.
-  E. Burnaev, M. Panov, and A. Zaytsev, Regression on the basis of nonstationary gaussian processes with bayesian regularization, Journal of communications technology and electronics, 61 (2016), pp. 661–671.
-  Y.-L. Chen, C.-T. C. Hsu, and H.-Y. M. Liao, Simultaneous tensor decomposition and completion using factor priors, IEEE Transactions on Pattern Analysis & Machine Intelligence, (2013), p. 1.
-  S. Gandy, B. Recht, and I. Yamada, Tensor completion and low-n-rank tensor recovery via convex optimization, Inverse Problems, 27 (2011), p. 025010.
-  L. Grasedyck, M. Kluge, and S. Krämer, Alternating directions fitting (adf) of hierarchical low rank tensors, Preprint, 149 (2013).
-  L. Grasedyck, M. Kluge, and S. Krämer, Alternating least squares tensor completion in the tt-format, arXiv preprint arXiv:1509.00311, (2015).
-  C.-Y. Ko, K. Batselier, W. Yu, and N. Wong, Fast and accurate tensor completion with total variation regularized tensor trains, arXiv preprint arXiv:1804.06128, (2018).
-  D. Kressner, M. Steinlechner, and B. Vandereycken, Low-rank tensor completion by riemannian optimization, BIT Numerical Mathematics, 54 (2014), pp. 447–468.
-  J. Liu, P. Musialski, P. Wonka, and J. Ye, Tensor completion for estimating missing values in visual data, IEEE transactions on pattern analysis and machine intelligence, 35 (2013), pp. 208–220.
-  M. Munkhoeva, Y. Kapushev, E. Burnaev, and I. Oseledets, Quadrature-based features for kernel approximation, in Advances in Neural Information Processing Systems, 2018, pp. 9147–9156.
-  I. Oseledets and E. Tyrtyshnikov, Tt-cross approximation for multidimensional arrays, Linear Algebra and its Applications, 432 (2010), pp. 70–88.
-  I. V. Oseledets, Tensor-train decomposition, SIAM Journal on Scientific Computing, 33 (2011), pp. 2295–2317.
C. E. Rasmussen, Gaussian processes in machine learning
, in Advanced lectures on machine learning, Springer, 2004, pp. 63–71.
-  B. Recht, M. Fazel, and P. A. Parrilo, Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization, SIAM review, 52 (2010), pp. 471–501.
-  M. Steinlechner, Riemannian optimization for high-dimensional tensor completion, SIAM Journal on Scientific Computing, 38 (2016), pp. S461–S484.
Low-rank tensor methods for linear systems and eigenvalue problems, PhD thesis, ETH Zurich, 2012.
-  E. Tyrtyshnikov, Incomplete cross approximation in the mosaic-skeleton method, Computing, 64 (2000), pp. 367–380.
-  T. Yokota, Q. Zhao, and A. Cichocki, Smooth parafac decomposition for tensor completion, IEEE Transactions on Signal Processing, 64 (2016), pp. 5423–5436.
-  L. Yuan, Q. Zhao, and J. Cao, Completion of high order tensor data with missing entries via tensor-train decomposition, in International Conference on Neural Information Processing, Springer, 2017, pp. 222–229.
-  Z. Zhang, T.-W. Weng, and L. Daniel, Big-data tensor recovery for high-dimensional uncertainty quantification of process variations, IEEE Transactions on Components, Packaging and Manufacturing Technology, 7 (2017), pp. 687–697.
-  L. Y. Zhao, J. Cao, et al., High-dimension tensor completion via gradient-based optimization under tensor-train format, arXiv preprint arXiv:1804.01983, (2018).
-  Q. Zhao, L. Zhang, and A. Cichocki, Bayesian cp factorization of incomplete tensors with automatic rank determination, IEEE transactions on pattern analysis and machine intelligence, 37 (2015), pp. 1751–1763.