Tensor Completion with Side Information: A Riemannian Manifold Approach

1 Introduction

Low Rank Tensor Completion (LRTC) problem, which aims to recover a tensor from its linear measurements, arises naturally in many artificial intelligence applications. In hyperspectral image inpainting, LRTC is applied to interpolate the unknown pixels based on the partial observation [Xu et al., 2015]. In recommendation tasks, LRTChelps users find interesting items [Liu et al., 2015]. In computational phenotyping, one adopts LRTC to discovery phenotypes in heterogeneous electronic health records [Wang et al., 2015].

Euclidean Models: LRTC can be formulated by a variety of optimization models over the Euclidean space. Amongst them, convex models that encapsulate LRTC as a regression problem penalized by a tensor nuclear norm are the most popular and well-understood [Romera-Paredes and Pontil, 2013; Zhang et al., 2014]. Though most of them have sound theoretical guarantees [Zhang and Aeron, 2016; Chen et al., 2013; Yuan and Zhang, 2015], in general, their solvers are ill-suited for large tensors because these procedures involve Singular Value Decomposition (SVD) of huge matrices per iteration [Liu et al., 2013]. Another class of Euclidean models is formulated as the decomposition problem that factorizes a low rank tensor into small factors [Jain and Oh, 2014;

∗Corresponding author

Filipovi´c and Juki´c, 2015; Xu et al., 2015]. Many solvers for such decomposition based model have been proposed, and low per-iteration computational cost is illustrated [Beutel et al., 2014; Liu et al., 2014; Smith et al., 2016].

Riemannian Models: LRTC can also be modeled by optimization constrained on Riemannian manifolds [Kressner et al., 2014; Kasai and Mishra, 2016], which is easily handled by many manifold based solvers [Absil et al., 2009]. Empirical comparison has shown that Riemannian solvers use significantly less CPU time to recover the underlying tensor in contrast to the Euclidean solvers [Kasai and Mishra, 2016]. The main reason resides in that they avoid SVD of huge matrices by explicitly exploiting the geometrical structure of LRTC, which makes them more suitable for massive problem.

Of all the Riemanian models, two search spaces, fix multilinear rank manifold [Kressner et al., 2014] and Tucker manifold [Kasai and Mishra, 2016], are usually employed. The former is a sub-manifold of Euclidean space, and the latter is a quotient manifold induced by the Tucker decomposition. Generally, quotient manifold based solvers have higher convergence rates because it is usually easier to design a pre-conditioner for them [Kasai and Mishra, 2016; Mishra and Sepulchre, 2016].

Side Information: In the Euclidean models of LRTC, side information has been proved to be helpful in improving the accuracy [Narita et al., 2011; Acar et al., 2011; Beutel et al., 2014]. One common form of the side information is the feature matrix, which measures the statistical properties of tensor modes [Kolda and Bader, 2009]. For example, in Netflix tasks, feature matrix can be built from the demography of users [Bell and Koren, 2007]. Another form is the similarity matrix, which quantifies the resemblance between two entities of a tensor mode. For instance, the social network generates the similarity matrix by utilizing the correspondence between users [Rai et al., 2015]. In practice, these two matrices can be transformed to each other, and we only consider the feature matrix case.

However, as far as we know, side information has not been incorporated in any Riemannian model for LRTC. The first difficulty lies in the model design. Fusing the side information into the Riemannian model inevitably compromises the integrity of the low rank tensor due to the compactness of the manifold. The second difficulty results from the solver design. Incorporating the side information may aggravate the ill-conditioning of LRTC problem and degenerates the convergence significantly.

Contributions: To address these difficulties, a novel Riemannian LRTC method is proposed from the perspective of both model and solver designs. By exploring the relation between the subspace spanned by the tensor fibers and the column space of the feature matrix, we explicitly integrate the side information in a compact way. Meanwhile, a first order solver is devised under the manifold optimization framework. To ease the ill-conditioning, we design a novel metric based on an approximated Hessian of the cost function. The metric implicitly induce an adaptive preconditioner for our solver. Empirical studies illustrate that our method achieves much more accurate solutions within comparable processing time than the state-of-the-art.

2 Notations and Preliminaries

In this paper, we only focus on the 3rd order tensor, but generalizing our method to higher order is straight forward. We use the notation X ∈ Rn×m to denote a matrix, and the notation X ∈ Rn1×···×nd to denote a d-th order tensor. We also denote by X(i1, · · · , id) the element in position (i1, · · · , id) of X. For many cases, we use abbreviation like {Oi}3i=1 to denote the sequence O1,O2,O3.

Mode-k Fiber and matricization: A fiber of a tensor is obtained by varying one index while fixing the others, i.e. X(i1, · · · , ik-1, :, ik+1, · · · , id) is the mode-k fiber of a d-th order tensor X. A mode-k matricization X(k) ∈ Rnk×(n1···nk-1nk+1···nd) of a tensor X is obtained by arranging the mode-n fibers of X so that each of them is a column of X(k) [Kolda and Bader, 2009] The mode-k product of tensor X and matrix A is denoted by X ×k A, whose mode-k matricization can be expressed as (X ×k A)(k) = AX(k). Inner product and norm: The inner product of two tensors with the same size is defined by 〈X,Y〉 = ∑i1,··· ,id X(i1, · · · , id)Y(i1, · · · , id). The Frobenius norm of a tensor X is defined by ‖X‖F = √〈X,X〉.

Multi-linear rank and Tucker decomposition: The multilinear rank rankvec (X) of a tensor X ∈ Rn1×n2×n3 is defined as a vector (rank(X(1)), rank(X(2)), rank(X(3))). If rankvec (X) = (r1, r2, r3), tucker decomposition factorizes X into a small core tensor G ∈ Rr1×r2×r3 and three matrices Ui ∈ Rni×ri with orthogonal columns, that is X = G×3i=1 Ui. Note that, the tucker decomposition of a tensor is not unique. In fact, if X = G ×3i=1 Ui, we can easily obtain X =H×3i=1 Vi, withH= G×3i=1 O⊤i , Vi =UiOi, where Oi ∈ Rri×ri is any orthogonal matrix. Thus, we obtain the equivalent class

We denote [G, {Ui}3i=1 ] by [X], when X = G ×3i=1 Ui. Usually, [X] is called the Tucker representation of X, while X is call the tensor representation of [X]. We also use X to denote a specific decomposition of X , additionally X ∈ [X ].

2.1 Search Space of RiemannianModels

The Tucker manifold that we used in our Riemannian model is a quotient manifold induced by the Tucker decomposition.

In order to lay the ground for Tucker manifold, we first describe its counterpart, the fix multi-rank manifold, which will be helpful in understanding the whole derivation.

A fixed multi-linear rank manifold Fr consists of tensors with the same fixed multi-linear rank. Specifically

To define the Tucker manifold, we first define a total space

in which S(ri, ni) is the Stiefel manifold of ni × ri matrices with orthogonal columns. Then, we can depict the Tucker manifold of multi-linear rank r as follows.

The Tucker manifold is a quotient manifold of the total space (1). We use the abstract quotient manifold, rather than the concrete total space, as search space because the nonuniqueness of the Tucker decomposition is undesirable for optimization. Note that such non-uniqueness will introduce more local optima into the minimization. The relation of manifold Fr and Mr/ ∼ is characterized as follows.

Proposition 1. The quotient manifold Mr/ ∼ is diffeomorphic to the fix multi-linear rank manifold Fr, with diffeomorphism ρ(·) from Fr to Mr/ ∼ defined by ρ(X) = [G, {Ui}3i=1] where [G, {Ui}3i=1] is the tucker representation of X.

This proposition says that each tensor X ∈ Fr can be represented by a unique equivalent class [G, {Ui}3i=1 ] ∈ M/ ∼ and vice-versa.

2.2 Vanilla Riemannian Tensor Completion

The purest incarnation of Riemannian tensor completion model is the Riemannian model over the fix multi-linear rank manifold. Let R ∈ Rn1×n2×n3 be a partially observed tensor. Let Ω be the set which contains the indices of observed entries. The model can be expressed as:

with PΩ maps X to the sparsified tensor PΩ(X), where PΩ(X)(i1, i2, i3) = X(i1, i2, i3) if (i1, i2, i3) ∈ Ω, and PΩ(X)(i1, i2, i3) = 0 otherwise.

Another popular model, Tucker model, is based on the quotient manifold Mr/ ∼, which can be expressed as:

with ρ defined in Prop. 1.

Note that since the dawn of Riemannian framework for LRTC, a quandary exists: on one hand, sparse measurement limits the capacity of the solution; on the other hand, rich side information can not be incorporated into this framework. In many artificial intelligence applications, demands for high accuracy further exacerbates such dilemma.

3 RiemannianModel with Side Information

We focus on the case that the side information is encoded in feature matrices Pi ∈ Rni×ki . Suppose R ∈ Fr has tucker factors (G, {Ui}3i=1 ).Without loss of generality, we assume that ki ≥ ri and Pi has orthogonal columns.

In the ideal case, we assume that

Such relation means that the feature matrices contain all the information in the latent space of the underlying tensor. Equivalently, there exists a matrix Wi such that Ui = PiWi. However, in practice, due to the existence of noise, one can only expect such relation to hold approximately, i.e. Ui ≈ PiWi. Incorporating such relation to a tensor completion model via penalization, we have the following formulation

where L(G, {Ui}3i=1 ) = ‖PΩ(G ×3i=1 Ui -R)‖2F /2. Fixing G and Ui, with respect toWi, (6) has a close form solution

Since minx,y l(x, y) = minx l(x, y(x)) where y(x) = argminy l(x, y), one can substitute (7) into the above problem and obtain the following equivalence

Although the cost function is already smooth over the total space Mr, due to its invariance over the equivalent class [G, {Ui}3i=1 ], there can be infinite local optima, which is extremely undesirable. Indeed, if (G, {Ui}3i=1 ) is a local optimal of the objective, then so is every point in the infinite set [G, {Ui}3i=1 ]. One way to reduce the number of local optima is to mathematically treated the entire set [G, {Ui}3i=1 ] as a point. Consequently, we redefine the cost by ˜ f([G, {Ui}3i=1 ]) = f(G, {Ui}3i=1 ) and obtain the following Remainnian optimization problem over the quotient manifold Mr/ ∼:min

Remark 1. In Riemannian optimization literature, problem (8) is called the lifted representation of problem (9) over the total space [Absil et al., 2009]. This model is closely related to the Laplace regularization model [Narita et al., 2011]. Concretely, they share the same form:

The difference lies in that Li is a projection matrix in our case, while, in the Laplace regularization model, Li is a Laplacian matrix.

Remark 2. Since each [X] ∈ Mr/ ∼ has a unique tensor representation in X ∈ Fr, we show that the abstract model (9) can be represented as a concrete model over the manifold Fr. Specifically, the following Proposition interprets the proposed model as an optimization problem with a regularizer that encourages the mode-i space of the estimated tensor close to span(Pi).

Proposition 2. if [X] is a critical point of problem (9) then its tensor representation X is a critical point of the following problem.

where dist(·, ·) is the Chodal distance [Ye and Lim, 2014] between two subspaces. And vice versa.

4 Riemannian Conjugate Gradient Descent

We depict the optimization framework for quotient manifolds in Fig. 1. Under this framework, we solve the proposed problem (9) by Riemannian Conjugate Gradient descent (CG). With the details specified later, we list our CG solver for problem (9) in Alg. 1, where the CG direction is composed in the Polak-Ribiere+ manner with the momentum weight β(k) computed by Flecher-Reeves formula [Absil et al., 2009], and Tk(·) is the projector of horizontal space HX (k) . The convergence property of CGSI (Alg. 1) to a stationary point follows the general analysis of CG method [Sato and Iwai, 2015]. To represent CGSI in concrete tensor formulations, four items must be specified: the Riemannian metric 〈·, ·〉X , the Riemannian gradient gradf(X), the retraction RX (·), and the projector onto horizontal space TX .

4.1 Metric Tuning

Riemannian metric 〈·, ·〉X of Mr is an inner product defined over each tangent space TX Mr. A high-quality Riemannian solver for a quotient manifold should be equipped with

Algorithm 1 CGSI: a Riemannian CG method

a well-tuned metric, because (1) the metric determines the differential structure of the quotient manifold, and more importantly (2) it implicitly endows the solver with a preconditioner, which heavily affects the convergent rate [Mishra and Sepulchre, 2014; Mishra, 2014].

From the perspective of preconditioning, it seems that the best candidate is the Newton metric 〈η, ξ〉X = D2 f(X)[η, ξ]∀η, ξ ∈ TX Mr where D2 f(X) is the second order differential of the cost function. However, under such metric, computing the search direction involves solving a large system of linear equations, which precludes the Newton metric from the application to huge datasets. Therefore, we propose to use the following alternative:

where g(X) is a scaled approximation to the original cost function, and D2 b g(X) is the block approximation of its second derivative, specifically g(X) ≜ 12 ‖G ×31 Ui -R‖2F + ∑3i=1 αiN 2 trace(UTi (Ii -PiPTi )Ui) with N = n1n2n3.

Our metric is more scalable than Newton metric. The following Proposition indicates that the scale gradient induced by this metric can be computed with O(∑3i=1 nikiri + r3 i ) additional operations.

Proposition 3. Suppose that the cost function f(·) has Euclidean gradient ∇f(X) = (∇Gf, {∇Ui f}3i=1). Then its scaled gradient ˜∇f(X) under the metric (11) can be computed by:

The final proposition suggests that the proposed metric makes the representation of solvers in the total space possible.

Proposition 4. The quotient manifoldMr/ ∼ admits a structure of Riemannian quotient manifold, if Mr is endowed with the Riemannian metric defined in (11).

4.2 Other Optimization Related Items

Projectors: To derive the optimization related items, two orthogonal projectors, ΨX (·) and ΠX (·), are required. The former projects a vector onto the tangent space TX Mr, and the latter is a projector from the tangent space onto the horizontal space HX . The orthogonality of both projectors is measured by the metric (11). For lack of space, the mathematical derivation is deferred to a long version of this paper.

Riemannian Gradient: According to [Absil et al., 2009], the Riemannian gradient can be computed by projecting the scaled gradient onto tangent space, specifically

Retraction: We use the retraction defined by

where uf(·) extracts the orthogonal component from a matrix. Such retraction is proposed by [Kasai and Mishra, 2016]. In the long version of this paper, we give rigorous analysis to prove that the above retraction is compatible with the proposed metric.

5 Experiments

We validate the effectiveness of the proposed solver CGSI by comparing it with the state-of-the-art. The baseline can be partitioned into three classes. The first class contains Riemannian solvers including GeomCG [Kressner et al., 2014], FTC [Kasai and Mishra, 2016], and gHOI [Liu et al., 2016]. The second class consists of Euclidean solvers that take no account of the side information, including AltMin [RomeraParedes et al., 2013] and HalRTC [Liu et al., 2013]. The third class comprises of two methods that incorporate side information, including RUBIK [Wang et al., 2015] and TFAI

[Narita et al., 2011]. All the experiments are performed in Matlab on the same machine with 3.0 GHz Intel E5-2690 CPU and 128GB RAM.

5.1 Hyperspectral Image Inpainting

A hyperspectral image is a tensor whose the slices are photographs of the same scene under different wavelets. We adopt the dataset provided in [Foster et al., 2006] which contains images about eight different rural scenes taken under 33 various wavelets. To make all methods in our baseline applicable to the completion problem, we resize each hyperspectral images to a small dimension such that n1 = 306, n2 = 402, and n3 = 33. Empirically, we treat these graphs as tensors of rank r = (30, 30, 6). The observed pixels, or the training set, are sampled from the tensors uniformly at random. And the sample size is set to |Ω| = OS×p in which OS is so-called Over-Sampling ratio and p = ∑3i=1 (niri - r2 i ) + r1r2r3 is the number of free parameters in a size n tensor with rank r. In addition to the observed entries, the mode-1 feature matrix is constructed by extracting the top(r1 + 10) singular vectors from a matrix of size n1 × 10r1 whose columns are sampled from the mode-1 fibers of the hyperspectral graphs. The recovery accuracy is measured by Normalized Root mean Square Error (NRSE) [Kressner et al., 2014]. All the compared methods are terminated when the training NRSE is less than 0.003 or iterate more than 300 epochs. We report the NRSE and CPU time of the compared methods in Tab. 2. From the table, we can see that the proposed method has much higher accuracy than the other solvers in our baseline. The empirical results also indicate that the sparser the observed pixels are the higher CGSI's improvement is on the recovery accuracy. The visual results of the 27th slices of recovered hyperspectral images of scene 7 are illustrated in Fig. 2.

5.2 Recommender System

In recommendation tasks, two datasets are considered: MovieLens 10M (ML10M) and MovieLens 20M (ML20M). Both datasets contain the rating history of users for items at specific moments. For both datasets, we partition the samples into 731 slices in terms of time stamp. Those slices have the identical time intervals. Accordingly, the completion tasks for the two datasets are of sizes 71567 × 10681 × 731 and 138493×26744×731 respectively. In addition to the rating history, both datasets contain two extra files: one describes the genres of each movie, and the other contains tags of each movie. We construct a corpus that contains the text description of all movies from the genres descriptions and all the tags. The feature matrix is extracted from the above corpus by the latent semantic analysis (LSA) method. The processing is efficient since LSA is implemented via randomized SVD.

Various empirical studies are conducted to validate the performance of the proposed method. In the first scenario, we record the CPU time and the Root Mean Square Error (RMSE) outputted by the compared algorithms under different choices of multi-linear rank. In this scenario, for both datasets, 80% samples are chosen as training set, and the rest are left for testing. The results are listed in Tab. 3, which suggests that the proposed method outperforms all other solvers in terms of accuracy. For ML10M, our method uses significantly less CPU time than its competitors. In Fig. 3, we report another scenario, in which the percentage of training samples are varied from 10% to 70% and the rank parameter is fixed to (10, 10, 10). Experimental results in this figure indicate that our method has the lowest RMSE.

To show the impact of parameter α on the performance of our method, we depict the relation between RMSE and α in Fig. 4, where the rank parameter is set to (10, 10, 10), and percentage of training samples is set to 80%. From this Figure we can see that our method has higher accuracy than the vanilla Riemannian model's solver FTC for a wide range of parameter choices.

6 Conclusion

In this paper, we exploit the side information to improve the accuracy of Riemannian tensor completion. A novel Riemannian model is proposed. To solve the model efficiently, we design a new Riemannian metric that implicitly induce an adaptive preconditioner for the solving procedure. Then, we devise a Riemannian conjugate gradient descent method using the well-tuned metric. Empirical results show that our solver outperforms the state-of-the-art.

Acknowledgments

This work is partially supported by National Natural Science Foundation of China (Grant No: 61472347, 61672376, and 61672449).