SPARTan: Scalable PARAFAC2 for Large & Sparse Data

TL;DR

SPARTan is a scalable, efficient algorithm for large-scale, sparse PARAFAC2 tensor decompositions, enabling analysis of complex, temporally-evolving medical data with significant performance improvements.

Contribution

We introduce SPARTan, a novel, scalable method that significantly improves the efficiency and memory usage of PARAFAC2 tensor decomposition for large, sparse datasets.

Findings

22X performance improvement over previous methods

Successfully applied to large, real-world medical datasets

Identified clinically meaningful, temporally-evolving phenotypes

Abstract

In exploratory tensor mining, a common problem is how to analyze a set of variables across a set of subjects whose observations do not align naturally. For example, when modeling medical features across a set of patients, the number and duration of treatments may vary widely in time, meaning there is no meaningful way to align their clinical records across time points for analysis purposes. To handle such data, the state-of-the-art tensor model is the so-called PARAFAC2, which yields interpretable and robust output and can naturally handle sparse data. However, its main limitation up to now has been the lack of efficient algorithms that can handle large-scale datasets. In this work, we fill this gap by developing a scalable method to compute the PARAFAC2 decomposition of large and sparse datasets, called SPARTan. Our method exploits special structure within PARAFAC2, leading to a novel algorithmic reformulation that is both fast (in absolute time) and more memory-efficient than prior work. We evaluate SPARTan on both synthetic and real datasets, showing 22X performance gains over the best previous implementation and also handling larger problem instances for which the baseline fails. Furthermore, we are able to apply SPARTan to the mining of temporally-evolving phenotypes on data taken from real and medically complex pediatric patients. The clinical meaningfulness of the phenotypes identified in this process, as well as their temporal evolution over time for several patients, have been endorsed by clinical experts.

Figures

Tables

Table 1. Table 1. Running time comparison: Time in minutes of one iteration for increasingly larger datasets (63m to 500m) and fixed target rank (two cases considered: R = { 10 , 40 } 𝑅 10 40 R=\{10,40\} ). The mode sizes for the datasets constructed are: 1 1 1 Mil. subjects, 5 5 5 K variables and a maximum of 100 100 100 observations per subject. OoM (Out of Memory) denotes that the execution failed due to the excessive amount of memory requested. Experiments are conducted on a server with 1TB of RAM.

Target Rank	10				40
#nnz(Millions)	63	125	250	500	63	125	250	500
SPARTan	7.4	8.9	11.5	15.4	14	18.4	61	114
Sparse PARAFAC2	24.4	60.1	72.3	194.5	275.2	408.1	OoM	OoM

Table 2. Table 2. Notations used throughout the paper.

Symbol	Definition
$𝒳,𝐗,𝐱,x$	Tensor, matrix, vector, scalar
${𝐗}^{†}$	Moore-Penrose pseudoinverse
$𝐗(:,i)$	Spans the entire $i$-th column of $𝐗$ (same for tensors)
$𝐗(i,:)$	Spans the entire $i$-th row of $𝐗$ (same for tensors)
$diag\left(𝐱\right)$	Diagonal matrix with vector $𝐱$ on the diagonal
$diag\left(𝐗\right)$	Extract diagonal of matrix $𝐗$
${𝐗}_k$	shorthand for $𝐗(:,:,k)$ (k-th frontal slice of tensor $𝒳$)
$\left\{{𝐗}_k\right\}$	the collection of ${𝐗}_k$ matrices, for all valid $k$
${𝐗}_{\left(n\right)}$	mode-$n$ matricization of tensor $𝒳$
$\circ$	Outer product
$\otimes$	Kronecker product
$\odot$	Khatri-Rao product
$\ast$	Hadamard (element-wise) product

Table 3. Table 3. Summary statistics for the real datasets of our experiments. K 𝐾 K is the number of subjects, J 𝐽 J is the number of variables, I k subscript 𝐼 𝑘 I_{k} is the number of observations for the k 𝑘 k -th subject and #nnz corresponds to the total number of non-zeros.

Dataset	$K$	$J$	$max\left(I_k\right)$	#nnz
CHOA	464,900	1,328	166	12.3 Mil.
MovieLens	25,249	26,096	19	8.9 Mil.

Table 4. Table 4. Subset of phenotypes discovered by PARAFAC2 corresponding to Figure  8 . The title annotation for each phenotype is provided by the medical expert. The red color corresponds to diagnoses and the blue color corresponds to medications.

Cancer	Weight
Chemotherapy	0.35
Leukemias [39.]	0.27
Immunity disorders [57.]	0.23
HEPARIN AND RELATED PREPARATIONS	0.6
ANTIEMETIC/ANTIVERTIGO AGENTS	0.34
SODIUM/SALINE PREPARATIONS	0.32
TOPICAL LOCAL ANESTHETICS	0.19
ANTIHISTAMINES - 1ST GENERATION	0.16
Neurological System Disorders	Weight
Other nervous system symptoms and disorders	0.56
Rehabilitation care; fitting of prostheses; and adjustment of devices [254.]	0.5
Residual codes; unclassified; all E codes [259. and 260.]	0.46
Other connective tissue disease [211.]	0.33
Other and unspecified metabolic; nutritional; and endocrine disorders	0.18

Equations

$$\mathcal{X}\approx \sum _{r=1}^R{\mathbf{u}}_r\circ {\mathbf{v}}_r\circ {\mathbf{w}}_r$$

$${\mathbf{X}}_k\approx \mathbf{U}{\mathbf{S}}_k{\mathbf{V}}^T$$

$${\mathbf{X}}_k\approx {\mathbf{U}}_k{\mathbf{S}}_k{\mathbf{V}}^T$$

$${\mathbf{U}}_k={\mathbf{Q}}_{k}H$$

$$\underset{\{{\mathbf{U}}_k\},\{{\mathbf{S}}_k\},\mathbf{V}}{\min }\sum _{k=1}^K||{\mathbf{X}}_k-{\mathbf{U}}_k{\mathbf{S}}_k{\mathbf{V}}^T|{|}_F^2$$

$$\underset{{\mathbf{Q}}_k}{\min }||{\mathbf{X}}_k-{\mathbf{Q}}_k\mathbf{H}{\mathbf{S}}_k{\mathbf{V}}^T|{|}_F^2$$

$$\underset{\mathbf{H},\{{\mathbf{S}}_k\},\mathbf{V}}{\min }\sum _{k=1}^K||{\mathbf{Q}}_k^T{\mathbf{X}}_k-\mathbf{H}{\mathbf{S}}_k{\mathbf{V}}^T|{|}_F^2$$

$${\mathbf{M}}^{(1)}={\mathbf{Y}}_{(1)}\left(\mathbf{W}\odot \mathbf{V}\right),$$

$$\mathbf{W}\odot\mathbf{V}=\left[\begin{array}[]{c}\mathbf{T}^{(1)}\\ \mathbf{T}^{(2)}\\ \vdots\\ \mathbf{T}^{(K)}\end{array}\right]$$

$${\mathbf{M}}^{(1)}=\sum _{k=1}^K{\mathbf{Y}}_k{\mathbf{T}}^{(k)}$$

$${\mathbf{T}}^{(k)}(i,:)=\mathbf{V}(i,:)\ast \mathbf{W}(k,:),$$

$$[{\mathbf{Y}}_k{\mathbf{T}}^{(k)}{]}_{j,:}=$$

$$\stackrel{Eq.\text{eqref}eq:mttkrp{1}_1}{=}$$

$$\stackrel{(a)}{=}$$

$$\stackrel{(b)}{=}$$

$${\mathbf{M}}^{(2)}={\mathbf{Y}}_{(2)}\left(\mathbf{W}\odot \mathbf{H}\right)$$

$${\mathbf{M}}^{(2)}=\sum _{k=1}^K{\mathbf{Y}}_k^T{\mathbf{T}}^{(k)}$$

$$[{\mathbf{Y}}_k^T{\mathbf{T}}^{(k)}{]}_{j,:}=\left({\mathbf{Y}}_k(:,j{)}^T\mathbf{H}\right)\ast \mathbf{W}(k,:)$$

$${\mathbf{M}}^{(3)}={\mathbf{Y}}_{(3)}\left(\mathbf{V}\odot \mathbf{H}\right)$$

$$\mathbf{M}^{(3)}(:,r)=\left[\begin{array}[]{c}{\mathbf{H}(:,r)}^{T}~{}\mathbf{Y}_{1}~{}\mathbf{V}(:,r)\\ \vdots\\ {\mathbf{H}(:,r)}^{T}~{}\mathbf{Y}_{K}~{}\mathbf{V}(:,r)\end{array}\right]$$

$${\mathbf{M}}^{(3)}(k,r)=$$

$$=$$

$${\mathbf{M}}^{(3)}(k,:)=dot\left(\mathbf{H},{\mathbf{Y}}_k\mathbf{V}\right)$$

Full text

In this work, we fill this gap by developing a scalable method to compute the PARAFAC2 decomposition of large and sparse datasets, called SPARTan. Our method exploits special structure within PARAFAC2, leading to a novel algorithmic reformulation that is both fast (in absolute time) and more memory-efficient than prior work. We evaluate SPARTan on both synthetic and real datasets, showing $22\times$ performance gains over the best previous implementation and also handling larger problem instances for which the baseline fails. Furthermore, we are able to apply SPARTan to the mining of temporally-evolving phenotypes on data taken from real and medically complex pediatric patients. The clinical meaningfulness of the phenotypes identified in this process, as well as their temporal evolution over time for several patients, have been endorsed by clinical experts.

Sparse Tensor Factorization, PARAFAC2, Phenotyping, Unsupervised learning

††copyright: rightsretained††doi: 10.475/123_4††isbn: 123-4567-24-567/08/06††conference: ACM conference; August 2017; Halifax - Nova Scotia, Canada††journalyear: 2017††price: 15.00††ccs: Information systems Data mining

1. Introduction

This paper concerns tensor-based analysis and mining of multi-modal data where observations are difficult or impossible to align naturally along one of its modes. A concrete example of such data is electronic health records (EHR), our primary motivating application. An EHR dataset contains longitudinal patient information, represented as an event sequence of multiple modalities such as diagnoses, medications, procedures, and lab results. An important characteristic of such event sequences is that there is no simple way to align observations in time across patients. For instance, different patients may have varying length records between the first admission and the most recent hospital discharge; or, two patients whose records’ have the same length may still not have a sensible chronological alignment as disease stages and patient progress vary.

For tensor methods, such data poses a significant challenge. Consider the most popular tensor analysis method in data mining, the canonical polyadic (CP) decomposition (also known as PARAFAC or CANDECOMP) (Hitchcock, 1927; Harshman, 1970; Carroll and Chang, 1970). A dataset with three modes might be stored as an $I\times J\times K$ tensor $\mathcal{X}$, which CP then decomposes into a sum of multi-way outer (rank-one) products, $\mathcal{X}~{}\approx~{}\sum_{r=1}^{R}\mathbf{u}_{r}\circ\mathbf{v}_{r}\circ\mathbf{w}_{r}$, where $\mathbf{u}_{r},\mathbf{v}_{r},\mathbf{w}_{r}$ are column vectors of size $I,J,K$, respectively, that effectively represents latent data concepts. Its popularity owes to its intuitive output structure and uniqueness property that makes the model reliable to interpret (Kruskal, 1977; Sidiropoulos and Bro, 2000; Kolda and Bader, 2009; Papalexakis et al., 2016; Sidiropoulos et al., 2016), as well as the existence of scalable algorithms and software (Bader and Kolda, 2007; Chi and Kolda, 2012; Bader et al., 2015). However, to make, say, the irregular time points in EHR one of the input tensor modes would require finding some way to align time. This fact is an inherent limitation of applying the CP model: any preprocessing to aggregate across time may lose temporal patterns (Ho et al., 2014a, b; Wang et al., 2015), while more sophisticated temporal feature extraction methods typically need continuous and sufficiently long temporal measures to work (Sun et al., 2008). Other proposed methods specific to healthcare applications may give some good results (Zhou et al., 2014; Wang et al., 2013, 2014) but lack the uniqueness guarantee; thus, it becomes harder to reliably extract the actual latent concepts as an equivalent arbitrary rotation of them will provide the same fit. All of these weaknesses apply in the EHR scenario outlined above.

In fact, the type of data in the motivating example are quite general: consider that we have $K$ subjects, for which we record $J$ variables and we permit each $k$-th subject to have $I_{k}$ observations, which are not necessarily comparable among the different subjects. For this type of data, Harshman proposed the PARAFAC2 model (Harshman, 1972b). It is a more flexible version of CP: while CP applies the same factors across a collection of matrices, PARAFAC2 instead applies the same factor along one mode and allows the other factor matrix to vary (Kolda and Bader, 2009). At the same time, it preserves the desirable properties of CP, such as uniqueness (Harshman and Lundy, 1996; ten Berge and Kiers, 1996; Kiers et al., 1999; Stegeman and Lam, 2015). As shown in Figure 1, PARAFAC2 approximates each one of the input matrices as: $\mathbf{X}_{k}\approx\mathbf{U}_{k}~{}\mathbf{S}_{k}~{}\mathbf{V}^{T}$, where $\mathbf{U}_{k}$ is of size $I_{k}\times R,\mathbf{S}_{k}$ is a diagonal $R$-by-$R$, $\mathbf{V}$ is of size $J\times R$ and $R$ is the target rank of the decomposition.

Despite its applicability, the lack of efficient PARAFAC2 decomposition algorithms has been cited as a reason for its limited popularity (Bro, 1997; Chew et al., 2007). Overall, PARAFAC2 has been mostly used for dense data (e.g., (Kiers et al., 1999)) or sparse data with a small number of subjects (Chew et al., 2007). To our knowledge, no work has assessed PARAFAC2 for large-scale sparse data, as well as the challenges arising by doing so.

In this paper, we propose SPARTan (abbreviated from Scalable PARafac Two) to fill this gap, with a focus on achieving scalability on large and sparse datasets. Our methodological advance is a new algorithm for scaling up the core computational kernel arising in the PARAFAC2 fitting algorithm. SPARTan achieves the best of both worlds in terms of speed and memory efficiency: a) it is faster than a highly-optimized baseline in all cases considered for both real (Figures 5, 6, 7) and synthetic (Table 1) datasets, achieving up to $22\times$ performance gain; b) at the same time, SPARTan is more scalable, in that it can execute in reasonable time for large problem instances when the baseline fails due to excessive memory consumption (Table 1).

We summarize our contributions as:

• •

Scalable PARAFAC2 method: We propose SPARTan, a scalable algorithm fitting the PARAFAC2 model on large and sparse data.

• •

Evaluation on various datasets: We evaluate the scalability of our approach using datasets originating from two different application domains, namely a longitudinal EHR and a time-evolving movie ratings’ dataset, which is also publicly available. Additionally, we perform synthetic data experiments.

• •

Real-world case study: We performed a case study of applying SPARTan on temporal phenotyping over medically complex pediatric patients in collaboration with Children’s Healthcare of Atlanta (CHOA). The phenotypes and temporal trends discovered were endorsed by a clinical expert from CHOA.

To promote reproducibility, our code will be open-sourced. For reviewing purposes, we provide the corresponding link 111https://goo.gl/ovQkQr.

2. Background

Next we describe the necessary terminology and operations regarding tensors. Then, we provide an overview of the CP model and relevant fitting algorithm. In Table 2, we summarize the notations used throughout the paper.

2.1. Tensors and Tensor Operations

The order of a tensor denotes the number of its dimensions, also known as ways or modes (e.g., matrices are $2$-order tensors). A fiber is a vector extracted from a tensor by fixing all modes but one. For example, a matrix column is a mode-$1$ fiber. A slice is a matrix extracted from a tensor by fixing all modes but two. In particular, the $\mathbf{X}(:,:,k)$ slices of a third-order tensor $\mathcal{X}$ are called the frontal ones and we succinctly denote them as $\mathbf{X}_{k}$ (Kolda and Bader, 2009). Matricization, also called reshaping or unfolding, logically reorganizes tensors into other forms without changing the values themselves. The mode-$n$ matricization of a $N$-order tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\dots\times I_{N}}$ is denoted by $\mathbf{X}_{(n)}\in\mathbb{R}^{I_{n}\times I_{1}I_{2}\dots I_{n-1}I_{n+1}\dots I_{N}}$ and arranges the mode-$n$ fibers of the tensor as columns of the resulting matrix.

2.2. CP Decomposition

The CP decomposition (Hitchcock, 1927; Harshman, 1970; Carroll and Chang, 1970) of a third-order tensor $\mathcal{X}\in\mathbb{R}^{I\times J\times K}$ is its approximation by a sum of three-way outer products:

[TABLE]

where $\mathbf{u}_{r}\in\mathbb{R}^{I},\mathbf{v}_{r}\in\mathbb{R}^{J}$ and $\mathbf{w}_{r}\in\mathbb{R}^{K}$ are column vectors. If we assemble the column vectors $\mathbf{u}_{r},\mathbf{v}_{r},\mathbf{w}_{r}$ as: $\mathbf{U}=[\mathbf{u}_{1}~{}\mathbf{u}_{2}~{}\dots~{}\mathbf{u}_{R}]\in\mathbb{R}^{I\times R},\mathbf{V}=[\mathbf{v}_{1}~{}\mathbf{v}_{2}~{}\dots~{}\mathbf{v}_{R}]\in\mathbb{R}^{J\times R},\mathbf{W}=[\mathbf{w}_{1}~{}\mathbf{w}_{2}~{}\dots~{}\mathbf{w}_{R}]\in\mathbb{R}^{K\times R}$, then $\mathbf{U},\mathbf{V},\mathbf{W}$ are called the factor matrices. Interpretation of CP is very intuitive: we consider that the input tensor can be summarized as $R$ latent concepts. Then, for each $r$-th concept, the vectors $(\mathbf{u}_{r},\mathbf{v}_{r},\mathbf{w}_{r})$ are considered as soft-clustering membership indicators, for the corresponding I, J and K elements of each mode. An equivalent formulation of Relation (1) w.r.t. the frontal slices $\mathbf{X}_{k}$ of the input tensor $\mathcal{X}$ is (Bro, 1998):

[TABLE]

where $k=1,2,\dots,K$ and $\mathcal{S}\in\mathbb{R}^{R\times R\times K}$ is an auxiliary tensor. Each frontal slice $\mathbf{S}_{k}$ of $\mathcal{S}$ contains the row vector $\mathbf{W}(k,:)$ along its diagonal: $\mathbf{S}_{k}=diag(\mathbf{W}(k,:))$. Relation (2) provides another viewpoint of interpreting the CP model, through its correspondence to the Singular Value Decomposition (SVD): each slice $\mathbf{X}_{k}$ is decomposed to a set of factor matrices $\mathbf{U},\mathbf{V}$ (similar to the singular vectors) which are common for all the slices, and a diagonal middle matrix (similar to the singular values) which varies for each $k$-th slice. Note, however, that no orthogonality constraints are imposed on $\mathbf{U},\mathbf{V}$ of the CP model, as in the SVD (Chew et al., 2007).

Uniqueness. A fundamental property of CP is uniqueness (Kruskal, 1977; Sidiropoulos and Bro, 2000). The issue with non-uniqueness can be exemplified via matrix factorization as follows (Kolda and Bader, 2009; Papalexakis et al., 2016): If a matrix $\mathbf{X}$ is approximated by the product of $\mathbf{A}\mathbf{B}^{T}$, then it can also be approximated with the same error by $\mathbf{A}\mathbf{Q}\mathbf{Q}^{-1}\mathbf{B}^{T}=\mathbf{\tilde{A}}\mathbf{\tilde{B}^{T}}$, for any invertible $\mathbf{Q}$. Thus, we can easily construct two completely different sets of rank-one factors that sum to the original matrix. Inevitably, this hurts interpretability, since we cannot know whether our solution is an arbitrarily rotated version of the actual latent factors. In contrast to matrix factorization or Tucker decomposition (Kolda and Bader, 2009), Kruskal (Kruskal, 1977) proved that CP is unique, under the condition: $k_{\mathbf{U}}+k_{\mathbf{V}}+k_{\mathbf{W}}\geq 2R+2$, where $k_{\mathbf{U}}$ is the $k$-rank of $\mathbf{U}$, defined as the maximum value $k$ such that any $k$ columns are linearly independent. The only exception is related to elementary indeterminacies of scaling and permutation of the component vectors (Kolda and Bader, 2009; Sidiropoulos et al., 2016). In sum, the CP decomposition is pursuing the true underlying latent information of the input tensor and provides reliable interpretation for unsupervised approaches.

Fitting the CP model. Perhaps the most popular algorithm for fitting the CP model is the CP-Alternating Least Squares (CP-ALS) (Carroll and Chang, 1970; Harshman, 1970), listed in Algorithm 1. The main idea is to solve for one factor matrix at a time, by fixing the others. In that way, each subproblem is reduced to a linear least-squares problem. In case the input tensor contains non-negative values, a non-negative least-squares solver (e.g., (Bro and De Jong, 1997)) can be used instead of an unconstrained one, to further improve the factors’ interpretability (Bro, 1997).

Due to the ever increasing need for CP decompositions in data mining, the parallel CP-ALS for sparse tensors has been extensively studied in the recent literature for both single-node and distributed settings (e.g., (Bader and Kolda, 2007; Kang et al., 2012; Choi and Vishwanathan, 2014; Papalexakis et al., 2015; Smith et al., 2015; Cheng et al., 2016)). A pioneering work in addressing scalability issues for sparse tensors was provided by Bader and Kolda (Bader and Kolda, 2007) 222The contributions of (Bader and Kolda, 2007), among others, are summarized as the Tensor Toolbox (Bader et al., 2015), which is widely acclaimed as the state-of-the-art package for single-node sparse tensor operations and algorithms.. The authors identified and scaled up the algorithm’s bottleneck, which is the materialization of the Matricized-Tensor-Times-Khatri-Rao-Product (MTTKRP). For example, in Algorithm 1, the MTTKRP corresponds to the computation of $\mathbf{X}_{(1)}(\mathbf{W}\odot\mathbf{V})$ when solving for $\mathbf{U}$. For large and sparse tensors, a naive construction of the MTTKRP requires huge storage and computational cost and has to be avoided.

3. PARAFAC2 Overview & Challenges

3.1. Model

As we introduced in Section 1, the PARAFAC2 model (Harshman, 1972b) can successfully deal with an incomparable mode of each slice $\mathbf{X}_{k}$ (Kiers et al., 1999). It does so, by introducing a set of $\mathbf{U}_{k}$ matrices replacing the $\mathbf{U}$ matrix of the CP model in Relation (2). Thus, each slice $\mathbf{X}_{k}$ is decomposed as shown in Figure 1:

[TABLE]

where $k=1,\dots,K,\mathbf{U}_{k}\in\mathbb{R}^{I_{k}\times R},\mathbf{S}_{k}\in\mathbb{R}^{R\times R}$ is diagonal and $\mathbf{V}\in\mathbb{R}^{J\times R}$. To preserve uniqueness, Harshman (Harshman, 1972b) imposed the constraint that the cross product $\mathbf{U}_{k}^{T}\mathbf{U}_{k}$ is invariant regardless which subject $k$ is involved (Kolda and Bader, 2009; Chew et al., 2007). In that way, the CP model’s invariance of the factor $\mathbf{U}_{k}$ itself (or $\mathbf{U}$ given its invariance to k), is relaxed (Acar and Yener, 2009). For the above constraint to hold, each $\mathbf{U}_{k}$ factor is decomposed as:

[TABLE]

where $\mathbf{Q}_{k}$ is of size $I_{k}\times R$ and has orthonormal columns, and $\mathbf{H}$ is an $R\times R$ matrix, which does not vary by $k$ (Kolda and Bader, 2009). Then, the constraint that $\mathbf{U}_{k}^{T}\mathbf{U}_{k}$ is constant over $k$ is implicitly enforced, as follows: $\mathbf{U}_{k}^{T}\mathbf{U}_{k}=\mathbf{H}^{T}\mathbf{Q}_{k}^{T}\mathbf{Q}_{k}\mathbf{H}=\mathbf{H}^{T}\mathbf{H}=\mathbf{\Phi}$.

There have been several results regarding the uniqueness property of PARAFAC2 (Harshman and Lundy, 1996; ten Berge and Kiers, 1996; Kiers et al., 1999). The most relevant towards the large-scale data scenario (i.e., the number of $K$ subjects can easily reach the order of hundreds of thousands), is that PARAFAC2 is unique for $K\geq 4$ (Stegeman and Lam, 2015).

3.2. Classical Algorithm for PARAFAC2

Below, we overview the classical algorithm for fitting PARAFAC2, proposed by Kiers et al. (1999). Their original algorithm expects dense data as input. Its objective function is as follows:

[TABLE]

subject to: $\mathbf{U}_{k}=\mathbf{Q}_{k}\mathbf{H},~{}\mathbf{Q}_{k}^{T}\mathbf{Q}_{k}=\mathbf{I}$ and $\mathbf{S}_{k}$ to be diagonal. Algorithm 2 follows an Alternating Least Squares (ALS) approach, divided in two distinct steps: first (lines 3-6), the set of column-orthonormal matrices $\{\mathbf{Q}_{k}\}$ is computed by fixing $\mathbf{H},\mathbf{V},\{\mathbf{S}_{k}\}$. This step can be derived by examining the pursuit of each $\mathbf{Q}_{k}$ as an individual Orthogonal Procrustes Problem (Golub and Van Loan, 2013) of the form:

[TABLE]

subject to $\mathbf{Q}_{k}^{T}\mathbf{Q}_{k}=\mathbf{I}$. Given the SVD of $\mathbf{H}\mathbf{S}_{k}\mathbf{V}^{T}\mathbf{X}_{k}^{T}$ as $\mathbf{P}_{k}\mathbf{\Sigma}_{k}\mathbf{Z}_{k}$, the minimum of objective (5) over column-orthonormal $\mathbf{Q}_{k}$ is given by $\mathbf{Q}_{k}=\mathbf{Z}_{k}\mathbf{P}_{k}^{T}$ (Kiers et al., 1999; Golub and Van Loan, 2013).

Second, after solving for and fixing $\{\mathbf{Q}_{k}\}$, we find a solution for the rest of the factors as:

[TABLE]

where $\mathbf{S}_{k}$ is diagonal. Note the equivalence of the above objective with the CP “slice-wise” formulation of Relation (2). This equivalence implies that minimizing the objective (6) is achieved by executing the CP decomposition on a tensor $\mathcal{Y}\in\mathbb{R}^{R\times J\times K}$ with frontal slices $\mathbf{Y}_{k}=\mathbf{Q}_{k}^{T}\mathbf{X}_{k}$ (lines 7-10). In order to avoid executing all the costly CP iterations, Kiers et al. (Kiers et al., 1999) propose to run a single CP-ALS iteration, since this suffices to decrease the objective.

The PARAFAC2 model can be extended so that non-negative constraints are imposed on $\{\mathbf{S}_{k}\},\mathbf{V}$ factors (Bro, 1998). This is a property inherited by the CP-ALS iteration, where we can constrain the factors $\mathbf{V}$ and $\mathbf{W}$ to be non-negative, as discussed in Section 2.2. Note that constraining the $\{\mathbf{U}_{k}\}$ factors to be non-negative as well is not as simple, and a naive approach would violate the model properties (Bro et al., 1999).

3.3. Challenges of PARAFAC2 on sparse data

Next, we summarize a set of crucial observations regarding the computational challenges, when executing Algorithm 2 on large, sparse data:

Bottleneck of Algorithm 2. Regarding the $1$st step (lines 3-6), in practice for sparse $\mathbf{X}_{k}$, each one of the $K$ sub-problems scales as $\mathcal{O}(min(RI^{2},R^{2}I))$, due to the SVD involved (Trefethen and Bau III, 1997), where $I$ is an upper bound for $I_{k}$. Note that this computation can be trivially parallelized for all $K$ subjects. On the other hand, the $2$nd step (lines 7-10) is dominated by the MTTKRP computation (which as we discussed in Section 2 is the bottleneck of sparse CP-ALS). Thus, it scales as $3R~{}nnz(\mathcal{Y})$ (Smith et al., 2015), using state-of-the-art sparse tensor libraries for single-node (Bader and Kolda, 2007). Given that none of the input matrices $\mathbf{X}_{k}$ is completely zero, then: $3R~{}nnz(\mathcal{Y})\geq 3KR^{2}$. As a result, this step becomes the bottleneck of Algorithm 2 for large and sparse “irregular” tensors, since it cannot be parallelized w.r.t. the $K$ subjects, as trivially happens with the $1$st one.

Imbalance of mode sizes of $\mathcal{Y}$. The size of the intermediate tensor $\mathcal{Y}$ formed is $R\times J\times K$. For large-scale data, we expect that $R<<K,J$ since $R$ corresponds to the target rank of the overall PARAFAC2 decomposition. Note that this property of “size imbalance” may not hold for general large-scale sparse tensors, but it is crucial for optimizing PARAFAC2.

Structured sparsity of $\{\mathbf{X}_{k}\}$. First, we observe that even if the input slices $\{\mathbf{X}_{k}\in\mathbb{R}^{I_{k}\times J}\}$ are very sparse, all their $I_{k}$ rows will contain at least one non-zero element. If this is not the case, we can simply filter the zero ones, without affecting the result. However, this does not hold for the $J$ columns of each $\mathbf{X}_{k}$, which have to be aligned across all $K$ subjects. A direct consequence of that can be easily noticed in real datasets. For example, in the EHR “irregular” tensor of observations $\times$ variables $\times$ subjects, very few variables (e.g., diagnostic codes) are recorded for each subject. Driven by this observation, we are motivated to computationally exploit any column sparsity (i.e., cases where many columns will be completely zero) of each one of the input matrices $\mathbf{X}_{k}$.

4. The SPARTan approach

4.1. Overview

Motivated by the challenges presented in Section 3, we propose a specialized version of the Matricized-Tensor-Times-Khatri-Rao-Product (MTTKRP) kernel, specifically targeting the intermediate tensor $\mathcal{Y}\in\mathbb{R}^{R\times J\times K}$ formed within the PARAFAC2-ALS algorithm. We first provide an overview of the properties that our approach exhibits:

• •

It is fully parallelizable w.r.t. the $K$ subjects. This property is crucial towards scaling up for large-scale “irregular” tensors.

• •

It exploits the structured sparsity of the input frontal slices $\{\mathbf{X}_{k}\}$. This is possible due to the observation that the $k$-th frontal slice $\mathbf{Y}_{k}=\mathbf{Q}_{k}^{T}\mathbf{X}_{k}$ of $\mathcal{Y}$ follows precisely the column sparsity pattern of $\mathbf{X}_{k}$. In particular, if $c_{k}$ is the number of columns of $\mathbf{X}_{k}$ containing at least one non-zero element, then $\mathbf{Y}_{k}$ will contain $R~{}c_{k}$ non-zero elements located in the positions of the non-zero columns of $X_{k}$. Exploiting structured sparsity is indispensable towards minimizing intermediate data and computations to the absolutely necessary ones.

• •

As a by-product of the above, SPARTan* avoids unnecessary data re-organization* (tensor reshaping/permutations), since all operations are formulated w.r.t. the frontal slices $\mathbf{Y}_{k}$ of tensor $\mathcal{Y}$. In fact, our approach never forms the tensor $\mathcal{Y}$ explicitly and directly utilizes the available collection of matrices $\{\mathbf{Y}_{k}\}$ instead.

4.2. Methodology

In the following, we describe the design of our MTTKRP kernel for each one of the tensor modes. We use the notation $\mathbf{M}^{(i)}$ to denote the MTTKRP corresponding to the $i$-th tensor mode. Note that our factor matrices are: $\mathbf{H}\in\mathbb{R}^{R\times R},\mathbf{V}\in\mathbb{R}^{J\times R}$ and $\mathbf{W}\in\mathbb{R}^{K\times R}$ as in Line 10 of Algorithm 2.

Mode-1 MTTKRP. First, we re-visit the MTTKRP equation:

[TABLE]

where $\mathbf{M}^{(1)}\in\mathbb{R}^{R\times R},\mathbf{Y}_{(1)}\in\mathbb{R}^{R\times KJ}$. In order to attempt to parallelize the above computation w.r.t. the $K$ subjects, we define the matrix $\mathbf{T}^{(k)}\in\mathbb{R}^{J\times R}$ to denote the $k$-th vertical block of the Khatri-Rao Product $\mathbf{W}\odot\mathbf{V}\in\mathbb{R}^{KJ\times R}$:

[TABLE]

We then remark that $\mathbf{Y}_{(1)}$ (i.e., mode-$1$ matricization of $\mathcal{Y}$) consists of an horizontal concatenation of the tensor’s frontal slices $\mathbf{Y}_{k}$. Thus, we exploit the fact that the matrix multiplication in Equation (7) can be expressed as the sum of outer products or more generally, as a sum of block-by-block matrix multiplications:

[TABLE]

Through Equation (8), the computation can be easily parallelized over $K$ independent sub-problems and then sum the partial results. This directly utilizes the frontal slices $\mathbf{Y}_{k}$ without further tensor organization. However, it constructs the whole Khatri-Rao Product (in the form of blocks $\mathbf{T}^{(k)}$). In order to avoid that, we first state an expression for each $i$-th row of $\mathbf{T}^{(k)}$, which is a direct consequence of the Khatri-Rao Product definition:

[TABLE]

where $*$ stands for the Hadamard (element-wise) product. Then, we express the $j$-th row of each partial result of Equation (8) as follows:

[TABLE]

where $(a)$ stems from the associative property of the Hadamard product and the fact that $\mathbf{W}(k,:)$ is independent of the summation and $(b)$ from the calculation of matrix multiplication as a sum of outer-products (in particular, we encounter the sub-case of vector-matrix product).

Equation (4.2) suggests an efficient way to compute the partial results of Equation (8), which we illustrate in Figure 2. First, we compute the matrix product $\mathbf{Y}_{k}\mathbf{V}$ and for each row of the intermediate result of size $R\times R$, we compute the Hadamard product with $\mathbf{W}(k,:)$. Note that, as we discussed in Section 3, $\mathbf{Y}_{k}$ is expected to be column-sparse in practice, thus multiplying by $\mathbf{V}$ uses only those rows of $\mathbf{V}$ corresponding to the non-zero columns of $\mathbf{Y}_{k}$. Thus, we avoid the redundant and expensive computation of the full Khatri-Rao Product. Overall, the methodology described above enjoys all of the properties described in Section 4.1.

Mode-2 MTTKRP. The methodology followed for the Mode-$2$ case is similar to the one described for the $1$st case. We state the corresponding MTTKRP equation:

[TABLE]

where $\mathbf{M}^{(2)}\in\mathbb{R}^{J\times R},\mathbf{Y}_{(2)}\in\mathbb{R}^{J\times RK}$. The main remark is that $\mathbf{Y}_{(2)}$ consists of an horizontal concatenation of the transposed frontal slices $\{\mathbf{Y}_{k}\}$ of the intermediate tensor $\mathcal{Y}$. Thus, if we denote as $\mathbf{T}^{(k)}$ the $k$-th vertical block of the Khatri-Rao Product $\mathbf{W}\odot\mathbf{H}$, we can formulate the problem as:

[TABLE]

Given the above, it is easy to extend Equation (4.2) for this case, so as to compute a single row of each partial result of Equation (12):

[TABLE]

The corresponding operations are illustrated in Figure 3. A crucial remark is that we can focus on computing the relevant intermediate results only for the non-zero rows of $\mathbf{Y_{k}}^{T}$, since the rest of the rows of the result $\mathbf{Y_{k}}^{T}\mathbf{T}^{(k)}$ will be zero. In sum, we again avoid redundant computations of the full Khatri-Rao Product and fulfill all of the properties described in Section 4.1.

Mode-3 MTTKRP. First, we state the equation regarding the Mode-$3$ case:

[TABLE]

where $\mathbf{M}^{(3)}\in\mathbb{R}^{K\times R}$ and $\mathbf{Y}_{(3)}\in\mathbb{R}^{K\times JR}$. Note that in this case, we are pursuing the MTTKRP of the mode corresponding to the $K$ subjects. Thus, an entirely different approach than the Mode-$1$, Mode-$2$ cases is needed so that we construct efficient independent sub-problems for each one of them. In particular, we need to design each $k$-th subproblem so that it computes the $k$-th row of $\mathbf{M}^{(3)}$. In addition, we want to operate only on $\{\mathbf{Y}_{k}\}$ without forming and reshaping the tensor $\mathcal{Y}$, as well as to exploit the frontal slices’ sparsity. To tackle the challenges above, we leverage the fact that (Choi and Vishwanathan, 2014):

[TABLE]

Then, we remark that in order to retrieve a certain element of the matrix $\mathbf{M}^{(3)}$, we have:

[TABLE]

The last line above reflects the inner product between the corresponding $r$-th columns of $\mathbf{H}$ and $[\mathbf{Y}_{k}~{}\mathbf{V}]$, respectively. Thus, in order to retrieve a row $\mathbf{M}(k,:)$, we can simply operate as:

[TABLE]

where the $dot()$ function extracts the inner product of the corresponding columns of its two matrix arguments. We illustrate this operation in Figure 4. Since $\mathbf{H}$ is a small $R$-by-$R$ matrix (due to the tensor’s “size imbalance”), it is very efficient to delay any computations on $\mathbf{H}$ until the $R$-by-$R$ intermediate matrix is formed as a product of $\mathbf{Y}_{k}~{}\mathbf{V}$. Then, we simply take the column-wise inner products between those two $R$-by-$R$ matrices. In that way, all the desirable properties we mentioned in Section 4.1 are also fulfilled.

In Algorithm 3, we list the pseudocode corresponding to the methodology proposed. Note that in lines 8,16, we can accumulate over the partial results in parallel, since the summation is independent of the iteration order.

5. Experiments

5.1. Setup

Real Data Description. Table 3 provides summary statistics regarding the real datasets used.

The CHOA (Children Healthcare of Atlanta) dataset corresponds to EHRs of pediatric patients with at least $2$ hospital visits. For each patient, we utilize the diagnostic codes and medication categories from their records, as well as the provided age of the patient (in days) at the visit time. The available International Classification of Diseases (ICD9) (Slee, 1978) codes are summarized to Clinical Classification Software (CCS) (ccs, 2017) categories, which is a standard step in healthcare analysis improving interpretability and clinical meaningfulness. We aggregate the time mode by week and all the medical events over each week are considered a single observation. The resulting data are of 464,900 subjects by 1,328 features by maximum 166 observations with 12.3m non-zeros.

MovieLens 20M is another real dataset we used, which is publicly available 333https://grouplens.org/datasets/movielens/. We are motivated to use this dataset, because of the importance of the evolution of user preferences over time, as highlighted in recent literature (Lathia et al., 2010). For this dataset, we consider that each year of ratings corresponds to a certain observation; thus, for each user, we have a year-by-movie matrix to describe her rating activity. We consider only the users having at least $2$ years of ratings.

Implementation details. We used MatlabR2015b for our implementations, along with functionalities for sparse tensors from the Tensor Toolbox (Bader et al., 2015) and the Non-Negative Least Squares (NNLS) approach (Bro and De Jong, 1997) from the N-way Toolbox (Andersson and Bro, 2000) 444We also accredit the dense PARAFAC2 implementation by Rasmus Bro, from where we have adapted many functionalities.. In both the SPARTan and the baseline implementations, we adjust the CP-ALS iteration arising in the PARAFAC2-ALS, so that non-negative constraints are imposed on the $\{\mathbf{S}_{k}\},\mathbf{V}$ factors, as discussed in Section 3.2.

The baseline method corresponds to the standard fitting algorithm for the PARAFAC2 model (Kiers et al., 1999) adjusted for sparse tensors as in (Chew et al., 2007). We utilized the implementation from the most recent version of the Tensor Toolbox (Bader et al., 2015) regarding both the manipulation of sparse tensors, as well as the CP-ALS iteration arising in the PARAFAC2-ALS.

Parallelism. We exploit the capabilities of the Parallel Computing Toolbox of Matlab, by utilizing its parallel pool in both SPARTan and the baseline approach, whenever this is appropriate. Regarding the size of the parallel pool, the number of workers of all the experiments regarding a certain dataset is fixed. For the movie-rating dataset we used the default of $12$ workers. For the synthetic and the CHOA datasets, we increased the number of workers to $20$ because of the data size increase.

Hardware. We conducted our experiments on a server running Ubuntu 14.04 with 1TB of RAM and four Intel E5-4620 v4 CPU’s with a maximum clock frequency of 2.10GHz. Each one of the processors contains $10$ cores, and each one of the cores can exploit $2$ threads with hyper-threading enabled.

5.2. SPARTan is fast and memory-efficient

Synthetic Data. We assess the scalability of the approaches under comparison for sparse synthetic data. We considered a setup with $1,000,000$ subjects, $5,000$ variables and a maximum of $100$ observations for each subject. The number of observations $I_{k}$ for each subject is dependent on the number of rows of $\mathbf{X}_{k}$ containing non-zero elements; thus, $I_{k}$ increases with the dataset density. Indicatively, the mean number of observations $I_{k}$ for the sparsest dataset created ($\approx 63$ mil.) is $46.9$ and for the densest ($\approx 500$ mil.) dataset, the mean $I_{k}$ is $99.3$. We randomly construct the factors of a rank-$40$ (which is the maximum target rank used in our experiments) PARAFAC2 model. Based on this model, we construct the input slices $\{\mathbf{X}_{k}\}$, which we then sparsify uniformly at random, for each sparsity level. The density of the sparsification governs the number of non-zeros of the collection of input matrices.

We provide the results in Table 1. First, we remark that SPARTan is both more scalable and faster than the baseline. In particular, the baseline approach fails to execute in the two largest problem instances for target rank $R=40$, due to out of memory problems, during the creation of the intermediate sparse tensor $\mathcal{Y}$. Note that as we discussed in Section 4, SPARTan avoids the additional overhead of explicitly constructing a sparse tensor structure, since it only operates directly on the tensor’s frontal slices $\{\mathbf{Y}_{k}\}$. Regarding the baseline’s memory issue, since the density of $\mathcal{Y}$ may grow (e.g., $\approx 10\%$ in the densest case), we also attempted to store the intermediate tensor $\mathcal{Y}$ as a dense one. However, this also failed, since the memory requested for a dense tensor of size $40$-by-$5$K-by-$1$Mil. exceeded the available RAM of our system ($1$TB). Overall, it is clear that the baseline approach cannot fully exploit the input sparsity. On the contrary, SPARTan properly executes for all the problem instances considered in a reasonable amount of time. In particular, for $R=40$, SPARTan is up to $22\times$ faster than the baseline. Even for a lower target rank of $R=10$, SPARTan achieves up to $13\times$ faster computation.

Real Data. We evaluate the scalability of the proposed SPARTan approach against the baseline method for the real datasets as well. In Figures 5, 6, 7, we present the results of the corresponding experiments. First, we target the full datasets and vary the pursued target rank (Figure 5). Note that for both datasets considered, the time per iteration of the baseline approach increases dramatically as we increase the target rank. On the contrary, the time required by SPARTan increases only slightly. Overall, our approach achieves up to over an order of magnitude gain regarding the time required per epoch for both datasets.

We also evaluate the scalability of the methods under comparison as we vary the subjects and the variables considered. Since the CHOA dataset (Figure 6) contains more subjects than variables, we vary the number of subjects for this dataset for two fixed target ranks ($10,40$). In both cases, SPARTan scales better than the baseline. As concerns the MovieLens dataset (Figure 7), since it contains more variables than subjects, we examine the scalability w.r.t. increasing subsets of variables considered. In this case as well, we remark the favorable scalability properties of SPARTan, rendering it practical to use for large and sparse “irregular” tensors.

5.3. Phenotype discovery on CHOA EHR Data

Motivation. Next we demonstrate the usefulness of PARAFAC2 towards temporal phenotyping of EHRs. Phenotyping refers to the process of extracting meaningful patient clusters (i.e., phenotypes) out of raw, noisy Electronic Health Records (Richesson et al., 2016). An open challenge in phenotyping is to capture temporal trends or patterns regarding the evolution of those phenotypes for each patient over time. Below, we illustrate how SPARTan can be used to successfully tackle this challenge.

Model Interpretation: We propose the following model interpretation towards the target challenge:

• •

The common factor matrix $\mathbf{V}$ reflects the phenotypes’ definition and the non-zero values of each $r$-th column indicate the membership of the corresponding medical feature to the $r$-th phenotype.

• •

The diagonal $\mathbf{S}_{k}$ provides the importance membership indicators of the $k$-th subject to each one of the $R$ phenotypes/clusters. Thus, we can sort the $R$ phenotypes based on the values of vector $diag(\mathbf{S}_{k})$ and identify the most relevant phenotypes for the $k$-th subject.

• •

Each $\mathbf{U}_{k}$ factor matrix provides the temporal signature of each patient: each $r$-th column of $\mathbf{U}_{k}$ reflects the evolution of the expression of the $r$-th phenotype for all the $I_{k}$ weeks of her medical history. Note that since all $\mathbf{X}_{k},\mathbf{S}_{k},\mathbf{V}$ matrices are non-negative, we only consider the non-negative elements of the temporal signatures in our interpretation.

Temporal Phenotyping of Medically Complex Patients (MCPs) In order to illustrate the use of PARAFAC2 towards temporal phenotyping, we focus our analysis on a subset of pediatric patients from CHOA, which are classified by them as Medically Complex. These are the patients with high utilization, multiple specialty visits and high severity. Conceptually, those patients suffer from chronic and/or very severe conditions that are hard to treat. As a result, it becomes a very important challenge to accurately phenotype those patients, as well as provide a temporal signature for each one of them, which summarizes their phenotypes’ evolution.

The number of MCPs in the CHOA cohort is $8,044$, their diagnoses and medications sum up to $1,126$, and the mean number of weekly observations for those patients is $28$. We ran SPARTan for target rank $R=5$ and the $2$ phenotypes corresponding to our example patient in Figure 8 are provided in Table 4 (phenotypes’ definition matrix). The rest of the phenotypes were also verified by the medical expert, but are omitted due to space limits. We provide the full phenotypes’ list, as well as additional patient examples in the Supplementary Material 555https://goo.gl/q7L91D. In Figure 8, we provide part of the real EHR, as well as the temporal signature produced by SPARTan, for a certain medically complex patient. Regarding the EHR, we visualize the subset of diagnoses and medications for which the sum of occurrences for the whole patient history is above a certain threshold (e.g., 5 occurrences). This step ensures that the visualized EHR will only contain the conditions exhibiting some form of temporal evolution. For the patient examples considered, we identify the top-$2$ relevant phenotypes through the importance membership indicator matrix $\mathbf{S}_{k}$ as discussed above. For those top-$2$ phenotypes, we present the resulting temporal signature, from which we easily detect intricate temporal trends of the phenotypes involved. Those trends were confirmed by the clinical expert as valuable towards fully understanding the phenotypic behavior of the MCPs.

6. Discussion & Conclusions

PARAFAC2 has been the state-of-the-art model for mining “irregular” tensors, where the observations along one of its modes do not align naturally. However, it has been highly disregarded by practitioners, as compared to other tensor approaches. Bro (Bro, 1997) has summarized the reason for that as:

The PARAFAC2 model has not yet been used very extensively maybe because the implementations so far have been complicated and slow.

The methodology proposed in this paper renders this statement no longer true for large and sparse data. In particular, as tested over real and synthetic datasets, SPARTan is both fast and memory-efficient, achieving up to $22\times$ performance gains over the best previous implementation and also handling larger problem instances for which the baseline fails due to insufficient memory.

The key insight driving SPARTan’s scalability is the pursuit and exploitation of special structure in the data involved in intermediate computations; prior art did not do so, instead treating those computations as a black-box.

The capability to run PARAFAC2 at larger scales is, in our view, an important enabling technology. As shown in our evaluations on EHR data, the clinically meaningful phenotypes and temporal trends identified by PARAFAC2 reflect the ease of the model’s interpretation and its potential utility in other application domains.

Future directions include, but are not limited to: a) development of PARAFAC2 algorithms for alternative models of computation, such as distributed clusters (Kang et al., 2012), or supercomputing environments; b) extension of the methodology proposed for higher-order “irregular” tensors with more than one mismatched mode.

Finally, to enable reproducibility and promote further popularization of the PARAFAC2 modeling within the area of data mining, we make our implementations publicly available.