VeST: Very Sparse Tucker Factorization of Large-Scale Tensors

TL;DR

VeST is a novel tensor factorization method that produces highly sparse and accurate decompositions of large-scale, partially observable tensors, improving interpretability and scalability over existing approaches.

Contribution

We introduce VeST, a scalable tensor factorization technique that yields very sparse core tensors and factor matrices without sacrificing accuracy.

Findings

VeST achieves at least 2.2 times more sparsity than competitors.

VeST is at least 2.8 times more accurate than existing methods.

The method scales well with input size and observable entries.

Abstract

Given a large tensor, how can we decompose it to sparse core tensor and factor matrices such that it is easier to interpret the results? How can we do this without reducing the accuracy? Existing approaches either output dense results or give low accuracy. In this paper, we propose VeST, a tensor factorization method for partially observable data to output a very sparse core tensor and factor matrices. VeST performs initial decomposition, determines unimportant entries in the decomposition results, removes the unimportant entries, and carefully updates the remaining entries. To determine unimportant entries, we define and use entry-wise 'responsibility' for the decomposed results. The entries are updated iteratively in a coordinate descent manner in parallel for scalable computation. Extensive experiments show that our method VeST is at least 2.2 times more sparse and at least 2.8 times more accurate compared to competitors. Moreover, VeST is scalable in terms of input order, dimension, and the number of observable entries. Thanks to VeST, we successfully interpret the result of real-world tensor data based on the sparsity pattern of the resulting factor matrices.

Tables

Table 1. Table 1 : Table of symbols and definitions.

$𝓧$	input tensor $\left(\in {ℝ}^{I_1\times \dots \times I_N}\right)$	$𝓖$	core tensor $\left(\in {ℝ}^{J_1\times \dots \times J_N}\right)$
$N$	order of $𝓧$	$I_n,J_n$	dimensionality of the $n$th mode of $𝓧$ and $𝓖$, respectively
${𝐀}^{\left(n\right)}$	$n$th factor matrix $\left(\in {ℝ}^{I_n\times J_n}\right)$	$a_{i_n{j}_n}^{\left(n\right)}$	$(i_n,j_n)$th element of ${𝐀}^{\left(n\right)}$
$\Omega$	set of observable entries of $𝓧$	$\left|\Omega \right|$	number of observable entries of $𝓧$
${\Omega }_{i_n}^{\left(n\right)}$	set of observable entries whose $n$th mode index is $i_n$	$\lambda$	regularization parameter for core and factor matrices
${\Vert 𝓧\Vert }_F$	Frobenius norm of tensor $𝓧$	${\Vert 𝓧\Vert }_1$	sum of absolute values of tensor $𝓧$
$\alpha$	an entry $(i_1,\dots ,i_N)$ of input tensor $𝓧$	$\beta$	an element $(j_1,\dots ,j_N)$ of core tensor $𝓖$
${\alpha }_{i_n=i}$	an entry $\left(i_1,\dots ,i_n=i,\dots ,i_N\right)$ of input tensor $𝓧$	${\beta }_{j_n=j}$	an element $\left(j_1,\dots ,j_n=j,\dots ,j_N\right)$ of core tensor $𝓖$

Table 2. Table 2 : Summary of datasets and hyperparameters used.

Name	Order	Dimensionality	Ranks	$\left|\Omega \right|$	${\left|\Omega \right|}_{test}$
MovieLens	4	$138K\times 27K\times 21\times 24$	$6\times 6\times 2\times 2$	18M	2M
Yelp	3	$71K\times 16K\times 108$	$10\times 10\times 10$	301K	33K
AmazonFood	3	$74K\times 256K\times 143$	$9\times 9\times 14$	511K	57K
Yelp-s	3	$50\times 50\times 10$	$5\times 5\times 5$	235	32
AmazonFood-s	3	$50\times 50\times 10$	$5\times 5\times 5$	444	51
Synthetic	$3-10$	${10}^3-{10}^8$	$3\times \cdots \times 3$	${10}^3-{10}^7$	-

Equations

$$L_F(\mathcal{G},{\mathbf{A}}^{(1)},...,{\mathbf{A}}^{(N)})=\sum _{\forall \alpha \in \Omega }{\left({\mathcal{X}}_{\alpha }-\sum _{\forall \beta \in \mathcal{G}}{\mathcal{G}}_{\beta }\prod _{n=1}^N{a}_{i_n{j}_n}^{(n)}\right)}^2+\lambda (\Vert \mathcal{G}{\Vert }_F^2+\sum _{n=1}^N{\Vert {\mathbf{A}}^{(n)}\Vert }_F^2),$$

$$L_1(\mathcal{G},{\mathbf{A}}^{(1)},...,{\mathbf{A}}^{(N)})=\sum _{\forall \alpha \in \Omega }{\left({\mathcal{X}}_{\alpha }-\sum _{\forall \beta \in \mathcal{G}}{\mathcal{G}}_{\beta }\prod _{n=1}^N{a}_{i_n{j}_n}^{(n)}\right)}^2+\lambda (\Vert \mathcal{G}{\Vert }_1+\sum _{n=1}^N{\Vert {\mathbf{A}}^{(n)}\Vert }_1)$$

$$Resp(e)=\frac{RE(e)-RE}{RE},$$

$$RE=\sqrt{\sum _{\forall \alpha =(i_1,\cdots ,i_N)\in \Omega }{\left({\mathcal{X}}_{\alpha }-\sum _{\forall \beta =(j_1,\cdots ,j_N)\in \mathcal{G}}{\mathcal{G}}_{\beta }\prod _{n=1}^N{a}_{i_n{j}_n}^{(n)}\right)}^2}/||\mathcal{X}|{|}_F$$

$$(RE(\bm{\mathscr{G}}_{\gamma}))^{2}=\frac{\sum\limits_{\forall\alpha=(i_{1},\cdots,i_{N})\in\Omega}\big{(}\bm{\mathscr{X}}_{\alpha}-\sum\limits_{\forall\beta=(j_{1},\cdots,j_{N})\neq\gamma\in\bm{\mathscr{G}}}\bm{\mathscr{G}}_{\beta}\prod\limits_{n=1}^{N}a^{(n)}_{i_{n}j_{n}}\big{)}^{2}}{||\bm{\mathscr{X}}||_{F}^{2}}$$

$$\begin{split}(RE(a^{(n)}_{ij}))^{2}=&{RE}^{2}+\frac{\sum\limits_{\forall\alpha\in\Omega_{i_{n}}^{(n)}}\big{(}2\cdot(\bm{\mathscr{X}}_{\alpha}-B(\alpha))+B_{j_{n}=j}(\alpha)\big{)}\cdot(B_{j_{n}=j}(\alpha))}{||\bm{\mathscr{X}}||_{F}^{2}},\end{split}$$

$${\mathcal{X}}_{\alpha =(i_1,\cdots ,i_N)}\approx B(\alpha )=\sum _{\forall \beta =(j_1,\cdots ,j_N)\in \mathcal{G}}{\mathcal{G}}_{\beta }\prod _{n=1}^N{a}_{i_n{j}_n}^{(n)},$$

$$B_{j_n=j}(\alpha )=\sum _{\forall {\beta }_{j_n=j}}{\mathcal{G}}_{{\beta }_{j_n=j}}\prod _{n=1}^N{a}_{i_n{j}_n}^{(n)},B_{j_n\ne j}(\alpha )=\sum _{\forall {\beta }_{j_n\ne j}}{\mathcal{G}}_{{\beta }_{j_n\ne j}}\prod _{n=1}^N{a}_{i_n{j}_n}^{(n)}.$$

$$a_{i_{n}j_{n}}^{(n)}\longleftarrow\underset{a_{i_{n}j_{n}}^{(n)}}{\operatorname{arg}\,\operatorname{min}}\;{\mathcal{L}_{F}(\bm{\mathscr{G}},\mathbf{A}^{(1)},...,\mathbf{A}^{(N)})}=\frac{\big{(}\sum\limits_{\forall\alpha\in\Omega_{i_{n}}^{(n)}}\bm{\mathscr{X}}_{\alpha}\delta_{\alpha}^{(n)}(j_{n})\big{)}-\big{(}\sum\limits_{\forall t\neq j_{n}}\bm{v}_{i_{n}j_{n}}^{(n)}(t)\cdot a_{i_{n}t}^{(n)}\big{)}}{\bm{v}_{i_{n}j_{n}}^{(n)}(j_{n})+\lambda},$$

$${\mathbf{v}}_{i_n{j}_n}^{(n)}(j)=\sum _{\forall \alpha \in {\Omega }_{i_n}^{(n)}}{\delta }_{\alpha }^{(n)}(j){\delta }_{\alpha }^{(n)}(j_n),$$

$${\delta }_{\alpha }^{(n)}(j)=\sum _{\forall {\beta }_{j_n=j}\in \mathcal{G}}{\mathcal{G}}_{{\beta }_{j_n=j}}\prod _{k\ne n}{a}_{i_k{j}_k}^{(k)},$$

$${\mathcal{G}}_{\beta }⟵\frac{\sum _{\forall \alpha \in \Omega }({\mathcal{X}}_{\alpha }-\sum _{\forall \gamma \ne \beta }{\mathcal{G}}_{\gamma }\prod _{n=1}^N{a}_{i_n{j}_n}^{(n)})\cdot \prod _{n=1}^N{a}_{i_n{j}_n}^{(n)}}{\lambda +\sum _{\forall \alpha \in \Omega }(\prod _{n=1}^N{a}_{i_n{j}_n}^{(n)}{)}^2}$$

$$\underset{a_{i_n{j}_n}^{(n)}}{\mathrm{argmin}}{L}_1(\mathcal{G},{\mathbf{A}}^{(1)},...,{\mathbf{A}}^{(N)})=\{\begin{array}{ll}(\lambda -g_{fm})/d_{fm}& \text{if gfm>λ}\\ -(\lambda +g_{fm})/d_{fm}& \text{if gfm<−λ}\\ 0& \text{otherwise}\end{array}$$

$$\underset{{\mathcal{G}}_{\beta }}{\mathrm{argmin}}{L}_1(\mathcal{G},{\mathbf{A}}^{(1)},...,{\mathbf{A}}^{(N)})=\{\begin{array}{ll}(\lambda -g_c)/d_c& \text{if g>λ}\\ -(\lambda +g_c)/d_c& \text{if g<−λ}\\ 0& \text{otherwise}\end{array}$$

Full text

1 Introduction

How can we factorize a large tensor to sparse core tensor and factor matrices such that outputs are easier to interpret? How can we do this without sacrificing accuracy? A tensor is a powerful tool for representing multi-modal data. Tensor factorization outputs a core tensor and factor matrices which reveal the latent relation of the data. Tensor factorization can also be viewed as a tool for multi-linear regression problem where only the target values, i.e., values of input tensor entries, are known. In this perspective, the columns of factor matrices act as latent components, their values as latent feature values, and the cells of the core tensor as weights of the relations between the latent components [8]. Sparse tensor factorization aims to output sparse core tensor and factor matrices. As sparse linear regression model enhances its interpretability [19], sparse factor matrices and a core improve interpretability.

There are two widely used approaches for sparse factorization. The first approach adds an $L_{1}$ norm as sparse regularizer to the factorization objective function [14, 12]. However, the sparsity is sensitive to lambda values. The second approach removes elements with small values from the core tensor or the factor matrices [21, 2, 22]. However, removing such elements does not necessarily leads to small reconstruction error, and thus value-based pruning sacrifices accuracy.

Other application-specific approaches include utilizing domain-specific knowledge as sparsity constraints [11], constructing factor matrix from sparse input sampling [10, 13], using smoothing matrices [17, 7], and learning sparse dictionary for image data [18, 23, 5]. These methods are based on several strong assumptions. The first approach assumes that prior classification of each mode, e.g., gene sets in omics data, are known. The second approach assumes that input tensors are already sparse and interpretable, e.g., network data. The third and fourth approaches assume values in input tensor are smooth, and unsmoothing them does not affect the overall meaning of the data, e.g., image data. However, these strong assumptions limit their use in general tensor factorizations.

In this paper, we propose VeST (Very Sparse Tucker factorization), a scalable and accurate Tucker factorization method to generate sparse factors and a core tensor for large-scale partially observable input tensor. VeST outputs very sparse factorization results by carefully determining the importance of elements of factors and the core, and pruning unimportant ones. VeST guarantees that the sparsity non-decreases in the update process by carefully derived update rules. Often, increasing sparsity too much degrades accuracy. VeST gives an algorithm to automatically determine a reasonable sparsity which offers a good balance with regard to accuracy. The very sparse result of VeST helps interpreting the result of Tucker factorization and easily revealing the relations of dimensions in multilinear regression.

Our main contributions are as follows:

• •

Algorithm. We propose VeST, an efficient Tucker factorization method for partially observable data, which produces very sparse outputs for better interpretability without loss of accuracy. VeST also provides an algorithm to automatically determine a sparsity which gives a good trade-off with regard to accuracy.

• •

Theory. We carefully derive parallelizable coordinate decent update rules for core tensor and factor matrices, and prove their correctness. We also analyze the time and space complexity of VeST.

• •

Performance. VeST provides better sparsity, accuracy, and scalability compared to other methods (see Figure 1).

The codes and datasets used in this paper are available at http://github.com/leesael/VeST.

2 Preliminaries and Related Works

We introduce concepts of tensor and its operations, Tucker factorization, and the standard algorithm for Tucker. Table 1 lists the symbols used.

2.1 Tensor and its Operations

Tensor is multi-dimensional array that contains numbers. An ‘order’ or ‘mode’ is the number of tensor dimensions, where a $1^{st}$-order tensor represents a vector and a $2^{nd}$-order tensor represents a matrix. We denote vectors by boldface lowercase letters (e.g., a), matrices by boldface capital letters (e.g., A), and three or higher order tensors by boldface Euler script letters (e.g., $\mathscr{X}$). An entry of a $3^{rd}$-order tensor can be expressed with three indices. For example, the ($i_{1}$, $i_{2}$, $i_{3}$)th entry of a $3^{rd}$-order tensor $\bm{\mathscr{X}}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$ is denoted by $x_{i_{1}i_{2}i_{3}}$, where index $i_{n}$ spans from 1 to $I_{n}$.

The size of a tensor is often evaluated by the Frobenius norm. Given an $N$-order tensor $\bm{\mathscr{X}}$ $(\in\mathbb{R}^{I_{1}\times...\times I_{N}})$, the Frobenius norm of $\bm{\mathscr{X}}$ is $||\bm{\mathscr{X}}||_{F}=\sqrt{\sum_{\forall{\alpha}\in\bm{\mathscr{X}}}{\bm{\mathscr{X}}^{2}_{\alpha}}}$, where $\alpha\ =(i_{1},\cdots,i_{N})$ is an index to an entry of input tensor $\bm{\mathscr{X}}$. Tensor decomposition often involves matricization of a tensor, and product between a tensor and a matrix. The mode-n matricization of a tensor $\bm{\mathscr{X}}\in\mathbb{R}^{I_{1}\times\cdots\times I_{N}}$ is denoted as $\mathbf{X}_{(n)}$ and the mapping from an entry $(i_{1},\cdots,i_{N})$ of $\bm{\mathscr{X}}$ to an entry $(i_{n},j)$ of $\mathbf{X}_{(n)}$ is given by $j=1+\sum_{k=1,k\neq n}^{N}[(i_{k}-1)\prod_{m=1,m\neq n}^{k-1}I_{m}]$. Also, the n-mode product of a tensor $\bm{\mathscr{X}}\in\mathbb{R}^{I_{1}\times\cdots\times I_{N}}$ with a matrix $\mathbf{U}\in\mathbb{R}^{J\times I_{n}}$ is denoted by $\bm{\mathscr{X}}\times_{n}\mathbf{U}$ ($\in\mathbb{R}^{I_{1}\times\cdots\times I_{n-1}\times J\times I_{n+1}\times\cdots\times I_{N}}$). Entry-wise, we have ${(\bm{\mathscr{X}}\times_{n}\mathbf{U})}_{i_{1}\cdots i_{n-1}ji_{n+1}\cdots i_{N}}=$ $\sum_{i=1}^{I_{n}}(\bm{\mathscr{X}}_{\alpha_{i_{n}=i}}{u}_{ji})$.

2.2 Tucker Factorization

Our proposed method VeST is built on top of Tucker factorization, one of the most popular tensor factorization methods. Given an $N^{th}$-order tensor $\bm{\mathscr{X}}$ $(\in\mathbb{R}^{I_{1}\times\cdots\times I_{N}})$, Tucker factorization approximates $\bm{\mathscr{X}}$ by a core tensor $\bm{\mathscr{G}}$ $(\in\mathbb{R}^{J_{1}\times\cdots\times J_{N}})$ and factor matrices $\{\mathbf{A}^{(n)}\in\mathbb{R}^{I_{n}\times J_{n}}|n=1\cdots N\}$ by minimizing the full reconstruction error: $\min_{\bm{\mathscr{G}},\mathbf{A}^{(1)},...,\mathbf{A}^{(N)}}\|\bm{\mathscr{X}}-\bm{\mathscr{G}}\times_{1}\mathbf{A}^{(1)}\cdots\times_{N}\mathbf{A}^{(N)}\|_{F}$.

Figure 2 illustrates a Tucker factorization result for a $3^{rd}$-order tensor. Typically, a core tensor $\bm{\mathscr{G}}$ is assumed to be smaller and denser than the input tensor $\bm{\mathscr{X}}$. Each factor matrix $\mathbf{A}^{(n)}$ represents the latent features of the object related to the $n$th mode of $\bm{\mathscr{X}}$, and each element of a core tensor $\bm{\mathscr{G}}$ indicates the weights of the relations composed of columns of factor matrices.

However, in real-world, data are often incomplete with some missing entries. To accommodate for the missing data, a partially observable Tucker factorization is needed. Given a tensor $\bm{\mathscr{X}}$ $(\in\mathbb{R}^{I_{1}\times...\times I_{N}})$ with observable entries $\Omega$, the goal of partially observable Tucker factorization of $\bm{\mathscr{X}}$ is to find factor matrices $\mathbf{A}^{(n)}$ $(\in\mathbb{R}^{I_{n}\times J_{n}},n=1,\cdots,N)$ and a core tensor $\bm{\mathscr{G}}~{}(\in\mathbb{R}^{J_{1}\times...\times J_{N}})$, which minimize the following loss:

[TABLE]

where $\alpha$ is an observable entry $(i_{1},...,i_{N})$ of input tensor $\bm{\mathscr{X}}$, $\beta$ is an element $(j_{1},...,j_{N})$ of core tensor $\bm{\mathscr{G}}$, and $\lambda>0$ is a regularization parameter. Note that the reconstruction error in Eq. (1) depends only on the observable entries of $\bm{\mathscr{X}}$, and $L_{F}$ regularization is used in Eq. (1) to prevent overfitting.

Tucker factorization often results in dense core and factor matrices. One of the approaches for sparsifying results is by including a sparsity constraint in the form of $L_{1}$ norm, a.k.a., Lasso, into the objective function. Given a tensor $\bm{\mathscr{X}}$ with observable entries $\Omega$, the goal of partially observable Tucker factorization via sparse regularizer of $\bm{\mathscr{X}}$ is to find factor matrices and a core tensor that minimize the following loss:

[TABLE]

which changed the $L_{2}$ regularization term of Eq. (1) to $L_{1}$. Again the reconstruction error in Eq. (2) depends only on the observable entries of $\bm{\mathscr{X}}$, and $L_{1}$ regularization is used to enforce sparsity. Another approach for sparsifying results is by pruning. That is, a partially observable Tucker factorization via minimal element value pruning of $\bm{\mathscr{X}}$ is to optimize on either Eq. (1) or Eq. (2), and sets the smallest $s$ ratio of the elements to zero in the core and factor matrices.

Evaluation of tensor decomposition and the prediction of the missing entry values (a.k.a., tensor completion) involves reconstruction. Given core tensor $\bm{\mathscr{G}}$ and factor matrices $\mathbf{A}^{(n)}$, the reconstruction of the original tensor $\bm{\mathscr{X}}$ is defined as $\bm{\mathscr{X}}\approx\bm{\mathscr{G}}\times_{1}\mathbf{A}^{(1)}\cdots\times_{N}\mathbf{A}^{(N)}$.

2.3 Tucker ALS Algorithm

A widely used technique for minimizing the loss functions Eq. (1) and Eq. (2) in a standard tensor factorization is alternating least squares (ALS) [8], which updates a factor matrix or a core tensor while keeping all others fixed.

Algorithm 1 describes a vanilla Tucker factorization algorithm based on ALS, which is called the higher-order orthogonal iteration (HOOI) (see [8] for details) that works on fully observable tensor. Notice that Algorithm 1 assumes missing entries of $\bm{\mathscr{X}}$ as zeros during the update process (lines 4-5). However, setting missing values to zero enforces Tucker ALS to factorize the original tensor such that missing values becomes zero when reconstructed. Note that the missing values are often nonzero values that are unknown. Thus, setting missing values to zero inserts false information into the factorization which results in higher reconstruction error as well as higher generalization error. Moreover, Algorithm 1 computes SVD (singular vector decomposition) given $\bm{\mathscr{Y}}_{(n)}$, which often results in dense matrices, thus tensor-ALS results in overall dense core tensor and factor matrices. Also, Algorithm 1 requires storing a full-dense matrix $\bm{\mathscr{Y}}_{(n)}$, and the amount of memory needed for storing $\bm{\mathscr{Y}}_{(n)}$ is $O(I_{n}\prod_{m\neq n}{J_{m}})$. The required memory grows rapidly when the order, the mode dimensionality, or the rank of a tensor increase, and ultimately causes intermediate data explosion [6].

In summary, the vanilla Tucker-ALS algorithm results in high generalization error in the presence of missing data, results in dense and thus hard-to-interpret core tensor and factor matrices, and cannot be applied to large data. Therefore, Algorithm 1 needs to be revised to focus only on observed entries, make sparse outputs, and be scaled for large-scale tensors at the same time.

3 Proposed Method

In this section, we propose VeST (Very Sparse Tucker factorization), a method for partially observed large scale tensor, that results in very sparse core tensor and factor matrices. Sparse results of tensor factorization increase interpretability and provide a scheme for better compression. To maximize sparsity without losing accuracy, VeST iteratively updates core tensor and factor matrices, and prunes unimportant elements from the core tensor and factor matrices. However, there are several challenges in designing an efficient update and pruning rules.

• •

Evaluating importance of elements. Vital elements of core tensor and factor matrices should not be pruned. How can we evaluate their importance?

• •

Automatically determining the sparsity. There is a trade-off relationship between sparsity and accuracy. How can we automatically determine an appropriate sparsity which gives a good balance with regards to accuracy?

• •

Updating factors while guaranteeing non-decreasing sparsity. The update process of factors and the core in the regular Tucker-ALS does not guarantee that the sparsity improves over the update process. How can we guarantee that update rules improve the sparsity?

We have the following main ideas to address the above challenges which we describe in detail in later subsections.

• •

Design responsibility indicator to evaluate contribution of each element on the accuracy (Section 3.2).

• •

Design auto-search algorithm VeST${}_{*}^{auto}$ to find a good sparsity that resides near the maximum sparsity just before the reconstruction error shoots up (Section 3.3).

• •

Design element-wise update rules to independently update each element of factor matrices and the core tensor. Element-wise update rules guarantee that the sparsity non-decreases by keeping pruned elements to zeros (Sections 3.4 and 3.5).

3.1 Overview

VeST is a scalable Tucker factorization method that results in very sparse core tensor and factor matrices for partially observable data (see Algorithm 2). First, VeST initializes all elements of the core tensor and factor matrices with random real values between 0 and 1 (line 1). Next, VeST iteratively updates the core tensor and factor matrices while pruning their elements (lines 3-7). In lines 3-4, VeST updates unpruned elements of the core tensor and factor matrices by element-wise update rules (Section 3.4), guaranteeing that the sparsity non-decreases. Then VeST prunes unimportant elements in the core tensor and factor matrices (lines 6-7). Importance of each element $e$ is evaluated by responsibility $Resp(e)$ which indicates how largely the element contributes to the accuracy (Section 3.2). $should\_prune()$ function determines when to stop pruning: if desired sparsity $s$ is achieved (in the manual version VeST${}_{*}^{man}$) or the reconstruction error shows a rapid increase (in the automatic version VeST${}_{*}^{auto}$). Motivated from simulated annealing, we gradually increase the pruning ratio as iterations proceed (line 7); this enables to explore larger search space in the beginning, while reducing the extent of the search to reduce to a minimum in the later iterations. The iterations proceed until the reconstruction error converges or the maximum iteration is reached. Finally, VeST standardizes all columns of factor matrices such that their norm is equal to one, and updates core tensor accordingly (lines 9-11). Specifically, $A^{(n)}$ is decomposed to $\mathbf{U^{(n)}}\in\mathbb{R}^{I_{n}\times J_{n}}$ and $\mathbf{B^{(n)}}$ where columns of $\mathbf{U^{(n)}}$ are unit vectors and $\mathbf{B^{(n)}}$ is a diagonal matrix whose $(i,i)$th element is the norm of $A^{(n)}$’s $i$th column. The core tensors are updated to maintain the same reconstruction error [9].

3.2 Evaluating Importance of Elements by Responsibility

VeST calculates the responsibility of each element which represents its contribution to the overall reconstruction accuracy over the observable elements of the input tensor to determine and prune unimportant elements. The intuition is that reconstruction error increases significantly when a vital element of the core tensor or factor matrices is set to zero, i.e., pruned. On the other hand, if the reconstruction error after pruning is similar or even smaller to that before pruning, the pruned element is insignificant. Formally, the responsibility is defined as follows.

Definition 1 (Responsibility)

Responsibility of an element $e$ in a factor matrix ($e=a^{(n)}_{ij}$) or core tensor ($e=\bm{\mathscr{G}}_{\gamma}$) is given by

[TABLE]

where

[TABLE]

is the normalized reconstruction error over the observable entries $\Omega$ of the original tensor $\bm{\mathscr{X}}$, and $RE(e)$ is residual reconstruction error defined when the element $e$ is set to zero (Eq. (5) and (6)). ∎

Definition 2 (Residual reconstruction error)

The residual reconstruction error $RE(\bm{\mathscr{G}}_{\gamma})$ for $(j_{1},...,j_{N})$th element $\gamma$ in core tensor $\bm{\mathscr{G}}$ is as follows:

[TABLE]

The residual reconstruction error $RE(a^{(n)}_{i,j})$ for an $(i,j)$th element $a^{(n)}_{i,j}$ in a factor matrix $\mathbf{A}^{(n)}$ is as follows:

[TABLE]

where $B(\alpha)$ is the entry-wise reconstruction defined as

[TABLE]

and $B_{j_{n}=j}(\alpha)$ and $B_{j_{n}\neq j}(\alpha)$ are the partial reconstruction functions defined as

[TABLE]

∎

The proof of correctness for the derivation of the Eq.(6) is provided in the supplementary material [16]. Note that both definitions are derived from the element-wise reformulation of the reconstruction error.

3.3 Pruning

After calculation of responsibility values, VeST prunes core tensor and factor matrices. Pruning is performed iteratively, each time after the core tensor and factor matrices are updated. The process of pruning an element consists of setting the value of the element to zero and marking the element as pruned in a marking table. The marked elements are excluded from the update step. To prune elements with low responsibility values, VeST sorts elements of the core tensor and each factor matrix, respectively, by the responsibility $Resp(e)$ in ascending order. Then, VeST prunes smallest $pr|\bm{\mathscr{G}}|$ elements from core tensor and smallest $pr|A^{(n)}|$ from each factor matrix, where $pr$ is the pruning rate of the current iteration. VeST starts with a small pruning rate $pr$ (INIT$\_$PR) and slowly increases $pr$ until maximum pruning rate (MAX$\_$PR) is reached. The default values of INIT$\_$PR and MAX$\_$PR are set to $0.01$ and $0.1$, respectively.

To determine when to stop pruning, VeST provides two different algorithms, the manual version VeST${}_{*}^{man}$ and the automatic version VeST${}_{*}^{auto}$; the subscript * denotes whether $L1$ or $L_{F}$ regularization is used.

• •

VeST${}_{*}^{man}$ takes a target sparsity $s$ as an input from the user and stops pruning when the total sparsity reaches $s$. That is, VeST${}_{*}^{man}$ enables users to decide on the lower bound of the final sparsity.

• •

VeST${}_{*}^{auto}$ determines the final sparsity automatically. VeST${}_{*}^{auto}$ does this by tracking changes in the reconstruction error and determines to stop pruning when an elbow point of the reconstruction error curve is reached. The elbow point is estimated as the point when the second derivative of the RE curve, estimated as $(RE_{t}+RE_{{t-2}}-2*RE_{{t-1}})/pr_{t}$ where $RE_{t}$ and $pr_{t}$ are RE and pruning rate at $t^{th}$ iteration, respectively, exceeds a small threshold (0.05 used).

3.4 Element-Wise Update Rules

VeST updates elements of the core tensor and factor matrices based on a coordinate decent approach in parallel. It enables VeST to update the core tensor and factor matrices without changing the value of the pruned elements. VeST checks the marking table which indicates whether elements have been pruned, and updates only the un-pruned elements. The update of an element is performed with observable tensor entries and fixed values of other elements in the factor matrices and the core tensor. The update rules for the core tensor and factor matrices are derived by setting the partial derivative of the loss function to zero and solving for each element. In previous works [15, 11], this approach has been shown and proven to converge faster with higher accuracy than existing approaches. Advantages of our update rules are that

1. accuracy is high and convergence is faster [15],
2. parallelization and selective updates are possible because all the elements are independently updated, and
3. the size of intermediate data is small, making the algorithm scalable.

3.4.1 Element-wise update rules with $L_{F}$ regularization

The update rule for an element $a_{i_{n}j_{n}}^{(n)}$ of factor matrix $A^{(n)}$ is derived by setting the partial derivative of loss function (Eq. (1)) with regard to $a_{i_{n}j_{n}}^{(n)}$ to zero.

Lemma 1 (Update rule for factor matrices with $L_{F}$ regularization)

[TABLE]

where $\bm{v}_{i_{n}j_{n}}^{(n)}$ is a length $J_{n}$ vector whose $j$th element is

[TABLE]

$\delta_{\alpha}^{(n)}$* is a length $J_{n}$ vector whose $j$th element is*

[TABLE]

$\Omega_{i_{n}}^{(n)}$* is the subset of $\Omega$ whose index of $n$th mode is $i_{n}$, and $\lambda$ is a regularization parameter. ∎*

The derivation of the core tensor update rule is similar to that of the factor matrix. The update rule for the $\beta^{th}$ element $\bm{\mathscr{G}}_{\beta}$ of the core tensor $\bm{\mathscr{G}}$ is given as follows.

Lemma 2 (Update rule for core tensor with $L_{F}$ regularization)

[TABLE]

∎

3.4.2 Element-wise update rules with $L_{1}$ regularization

For an element $a_{i_{n}j_{n}}^{(n)}$ of factor matrix $A^{(n)}$, the element-wise update rule with $L_{1}$ regularization is provided in the following Lemmas.

Lemma 3 (Update rule for factor matrix with $L_{1}$ regularization)

[TABLE]

where g_{fm}=2\big{(}\sum\limits_{\forall\alpha\in\Omega_{i_{n}}^{(n)}}\bm{\mathscr{X}}_{\alpha}\delta_{\alpha}^{(n)}(j_{n})\big{)}-\big{(}\sum\limits_{\forall t\neq j_{n}}\bm{v}_{i_{n}j_{n}}^{(n)}(t)\cdot a_{i_{n}t}^{(n)}\big{)}, $d_{fm}=2\bm{v}_{i_{n}j_{n}}^{(n)}(j_{n})$, and $\bm{v}_{i_{n}j_{n}}^{(n)}$, $\delta_{\alpha}^{(n)}$, $\Omega_{i_{n}}^{(n)}$, and $\lambda$ follow the same specification provided in Lemma 1. ∎

For an element $\bm{\mathscr{G}}_{\beta}$ of core tensor, the element-wise update rule with $L_{1}$ regularization is as follows:

Lemma 4 (Update rule for core tensor with $L_{1}$ regularization)

[TABLE]

where g_{c}=-2\sum\limits_{\forall\alpha\in\Omega}\big{(}\bm{\mathscr{X}}_{\alpha}-\sum\limits_{\forall\gamma\neq\beta}\bm{\mathscr{G}}_{\gamma}\prod\limits_{n=1}^{N}a_{i_{n}j_{n}}^{(n)}\big{)}\cdot\prod\limits_{n=1}^{N}a_{i_{n}j_{n}}^{(n)}, and d_{c}=2\sum\limits_{\forall\alpha\in\Omega}\big{(}\prod\limits_{n=1}^{N}a_{i_{n}j_{n}}^{(n)}\big{)}^{2}. ∎

The proofs of Lemmas 1 to 4 are provided in the supplementary material [16].

3.5 Parallel Update Algorithms

Responsibility calculation and factor matrices updates are performed in parallel. Algorithm 3 describes the pruning process where responsibility values of the core tensor and factor matrices are calculated in parallel for each observable entries of the input tensor. Note that the use of $B(\alpha)$ in line 5 enabled fast computing of $Resp(\bm{\mathscr{G}}_{\beta})$ in line 6; for a given $\beta$ in line 4, computing line 5 requires $O(|\Omega|)$ rather than $O(|\Omega||\bm{\mathscr{G}}|)$ since there is no need to compute $\sum_{\forall\beta\neq\gamma\in\bm{\mathscr{G}}}{\bm{\mathscr{G}}}_{\gamma}\prod_{n=1}^{N}a^{(n)}_{i_{n}j_{n}}$ in Eq. (5) from scratch.

The element-wise update of factor matrix $A^{(n)}$ is performed in parallel for each rows of factor matrices using either the $L_{F}$ or $L_{1}$ regularization (see Algorithm 4). Elements of the core tensor are dependent on each other and thus cannot be updated in parallel. However, considering that typical size ${|\bm{\mathscr{G}}|}$ of the core tensor is small, the core tensor updates are a minor burden in the computational process.

Lemma 5 (Complexity of VeST per iteration)

The time complexity per iteration of VeST is $O(N^{2}J|\bm{\mathscr{G}}||\Omega|/T+NIJ\log(IJ)+|\bm{\mathscr{G}}|\log|\bm{\mathscr{G}}|)$, and the memory complexity is $~{}O(TJ+J^{N}+NIJ)~{}$, where T is the number of threads, N is the order, I is the dimensionality of input tensor, and J is the dimensionality of the core tensor assuming the dimensionalities are equal for all modes. ∎

The full proof of Lemma 5 is in the supplementary material [16].

4 Experiments

We conduct experiments to answer the following questions.

Performance comparison (Section 4.2). How accurately and sparsely does VeST decompose a given tensor compared to other methods? 2. 2.

Sparsity and accuracy (Section 4.3). Does VeST successfully prune redundant information in the decomposition without hurting the accuracy? Does VeST${}_{*}^{auto}$ find reasonable sparsity and accuracy trade-off point? 3. 3.

Data scalability (Section 4.4). How scalable is VeST? 4. 4.

Interpretable discoveries (Section 4.5). How interpretable are the VeST results for discoveries on real-world tensors?

4.1 Experimental Settings

Datasets. We used three real-world datasets and synthetic datasets as summarized in Table 2. The real-world datasets are MovieLens111https://grouplens.org/datasets/movielens/, Yelp222http://www.yelp.com/dataset_challenge/, and AmazonFood333 http://snap.stanford.edu/data/web-FineFoods.html. MovieLens is a $4^{th}$ order tensor of movie ratings containing (user, movie, year, hour). Yelp is a $3^{rd}$ order tensor of business services rating data containing (user, business, year-month). AmazonFood is a $3^{rd}$ order tensor of food review scores from Amazon containing (product, user, year-month). To compare with other methods, we used subsets Yelp-s and AmazonFood-s of 3rd order tensors which are made denser than their originals. The density of Yelp-s and AmazonFood-s are 0.01 and 0.02, respectively. We also generated synthetic random tensors of various sizes and orders to test data scalability.

Environment. VeST was written in C++ with OpenMP [4] and Armadillo [20] libraries for parallelization. Methods L1 (Lasso) and Value Pruning were run on VeST framework with the difference just in the pruning approaches. We used the codes provided by the authors for TTP [21] (R) and Sparse CP [2] (Matlab). Tucker-ALS was performed via Tensor Toolbox for Matlab [3]. All experiments were done on a single machine equipped with an Intel Xeon E5-2630 v4 2.2GHz CPU (10 cores/20 threads) and 512GB memory. All reported measures are averages of five runs, unless otherwise stated.

Competitors. We compared VeST with the following methods.

1. $\bullet$

L1 (Lasso): A Tucker factorization method with lasso sparsity constraint implemented as VeST${}^{man}_{L_{1}}$ with sparsity $s=0$. 2. $\bullet$

Value Pruning: A Tucker factorization method with value pruning at the last step implemented as VeST${}_{*}^{man}$ with sparsity $s=0$ followed by value pruning with ratio 0.6 for $L_{F}$ and $L_{1}$ losses. 3. $\bullet$

TTP [21]: A tensor decomposition method that results in sparse components. 4. $\bullet$

Sparse CP [2]: CP decomposition method with lasso penalty. 5. $\bullet$

Tucker-ALS [8]: Conventional Tucker factorization method (HOOI). 6. $\bullet$

PTucker-Approx [15] : Tucker decomposition method for partially observable tensor with iterative element value pruning.

4.2 Performance Comparison

We compared the accuracy of VeST${}^{auto}_{L_{F}}$ and VeST${}^{auto}_{L_{1}}$ with those of the competitors on datasets Yelp-s and AmazonFood-s (Table 2). The comparison was performed on smaller and denser datasets of order three and not on original real-world datasets due to limitations of the TTR and Sparse CP.

We measured and compared normalized reconstruction errors (RE) over observable entries in input tensors. As shown in Fig. 1(a), VeST${}^{auto}_{L_{F}}$ and VeST${}^{auto}_{L_{1}}$ decomposed a given tensor with at least 2.8 times lower RE compared to other methods at a similar sparsity. Fig. 1(a) also shows that at a similar RE value, VeST${}_{L_{F}}^{*}$ and VeST${}_{L_{1}}^{*}$ output at least 2.2 times more sparse factor matrices and core tensor compared to other methods, where the sparsity is measured as the ratio of number of nonzero values in $\bm{\mathscr{G}}$ and $\mathbf{A}^{(n)}$ over $|\bm{\mathscr{G}}|+|\sum_{n=1}^{N}{\mathbf{A}^{(n)}}|$.

To answer how well VeST predicts missing entries compared to other methods, we measured REs of the reconstructed missing values (Test RE). After learning the factor matrices and the core tensor using 90% of the observed entries, we calculated the Test REs on the remaining $10\%$ of the observable entries. Fig 1(b) shows that VeST predicts missing entries at least 1.8 times more accurately compared to others, in addition to providing at least 1.7 times sparser results.

4.3 Sparsity and Accuracy

We tested the sparsity and accuracy of outputs of VeST on three full size real-world datasets: MovieLens, Yelp, and AmazonFood, in Figure 3. First, we investigated how the sparsity affects the normalized reconstruction error (RE). Note that the REs are not affected much until the sparsities are above 0.6; this shows that there is redundant information in the decomposition results, and VeST successfully finds and removes such redundancies to get compact outputs. Second, we investigated how VeST${}_{*}^{auto}$ automatically finds a desired sparsity which gives a good tradeoff with regard to accuracy. Note that in all the datasets, VeST${}_{*}^{auto}$ successfully finds sparsity points at the elbows of the RE curves, resulting in reasonable sparsity and RE trade-offs. Similar trend was observed for the sparsity and the test RE (see the supplementary material[16]).

4.4 Data Scalability

We evaluated the scalability of VeST by generating synthetic tensors varying the order, the dimensionality, and the number of observable entries, and measuring the running time (see Figure 4). For convenience, the dimensionality of each mode in the input tensor, as well as the dimensionality of core tensor, were set equal, i.e., $I_{1}=I_{2}=...=I_{N}$ and $J_{1}=J_{2}=...=J_{N}=3$, respectively.

Order. Data scalability on the order of input tensor is tested on synthetic tensors of varying orders from 3 to 10. For each input tensor, the dimensionality and the number of observable entries were fixed to $|\Omega|=10^{3}$. Figure 4 (a) shows that VeST scales quadratically with regard to order, as discussed in Section 3.5.

Dimensionality. Data scalability on the dimensionality was tested on input tensors of varying dimensionalities from $10^{3}$ to $10^{8}$ with each mode having equal dimension. The order was set to three, and $|\Omega|$ was set equal to the dimensionality. Figure 4 (b) shows that VeST has near-linear scalability in terms of the dimensionality.

Number of Observable Entries. Data scalability on the number of observable entries was tested on input tensors by varying the number of observable entries from $10^{3}$ to $10^{7}$. The order was set to three, and the dimensionality was fixed to $10^{3}$. Figure 4 (c) shows that VeST has near-linear scalability in term of the number of observable entries.

Effectiveness of Parallelization. We evaluated the parallelization scalability of VeST by increasing the number of threads from 1 to 20 and measuring $Time_{1}/Time_{T}$ where $Time_{T}$ is the running time per iteration using $T$ threads. Figure 4(d) shows near-linear scalability of VeST in terms of the number of threads used.

4.5 Discovery

We evaluated interpretability of VeST by investigating the factorization results of MovieLens dataset and visually showing that the sparse results enhance interpretability. It is difficult to analyze dense results generated by vanilla methods without post-processing. In contrast to existing methods, we can easily identify interesting factors generated by VeST based on sparsity of each row of a factor matrix.

Discovery of Greatest Movies. We found that a few rows of the movie-associated factor matrix are fully dense although the goal of VeST is to generate sparse results. Such rows corresponded to popular movies rated by diverse users. Figure 5 shows popular movies which have the largest number of non-zero and the sums of values. According to Empire magazine [1], 14 out of 20 movies we found were included in the 100 greatest movies. The remaining six movies, including ‘Sixth Sense’, ‘Kill Bill’, and ‘Fifth Element’, that were not in the 100 list were also very famous.

5 Conclusion

We proposed VeST, a very-sparse Tucker factorization method for sparse and partially observable tensors. By deriving the element-wise partial differential equations, determining the importance of elements by responsibilities, and parallel distribution of computational work, VeST successfully offers very sparse and accurate results that are applicable for large partially observable tensors. VeST generates at least 2.2 times more sparse results compared to other methods for partially observable tensors and at least 2.8 times accurate results compared to other sparse factorization methods. VeST also shows near linear scalability regarding tensor dimensionality, number of observable entries, and number of threads. Thanks to the increased sparsity that leads to improved interpretability by VeST, we were able to discover interesting patterns related to the greatest movies in the factor matrix of a real-world movie rating tensor data. Future works include better initialization for Tucker factorization, integration of prior knowledge, and effective visualization of tensor results.