Accelerating jackknife resampling for the Canonical Polyadic Decomposition

For vectors and matrices, we use bold lowercase and uppercase roman letters, respectively, e.g., $\mathbf{v}$ and $\mathbf{U}$. For tensors, we follow the notation in (Kolda and Bader, 2009); specifically, we use bold calligraphic fonts, e.g., $\boldsymbol{\mathcal{T}}$. The order (number of indices or modes) of a tensor is denoted by uppercase roman letters, e.g., $N$. For each mode $n\in\{1,2,\ldots,N\}$, a tensor $\boldsymbol{\mathcal{T}}$ can be unfolded (matricized) into a matrix, denoted by $\mathbf{T}_{(n)}$, where the columns are the mode-$n$ fibers of $\boldsymbol{\mathcal{T}}$, i.e., the vectors obtained by fixing all indices except for mode $n$. Sets are denoted by calligraphic fonts, e.g., $\mathcal{S}$. Given two matrices $\mathbf{A}$ and $\mathbf{B}$ with the same number of columns, the Khatri-Rao product, denoted by $\mathbf{A}\odot\mathbf{B}$, is the column-wise Kronecker product of $\mathbf{A}$ and $\mathbf{B}$. Finally, the unary operator $\oplus$, when applied to a matrix, denotes the scalar which is the sum of all matrix elements.

The standard alternating least-squares method for CP is shown in Algorithm 1 (CP-ALS). The input consists of a target tensor $\boldsymbol{\mathcal{T}}$. The output consists of a CP model represented by a sequence of factor matrices $\mathbf{U}_{1}$, $\ldots$ , $\mathbf{U}_{N}$. The algorithm repeatedly updates the factor matrices one by one in sequence until either of the following criteria are met: a) the fit of the model to the target tensor falls below a certain tolerance threshold, or b) a maximum number of iterations has been reached. To update a specific factor matrix $\mathbf{U}_{n}$, the gradient of the least-squares objective function with respect to that factor matrix is set to zero and the resulting linear least-squares problem is solved directly from the normal equations. This entails computing the Matricized Tensor Times Khatri-Rao Product (MTTKRP) (line 1), which is the product between the mode-$n$ unfolding $\mathbf{T}_{(n)}$ and the Khatri-Rao Product (KRP) of all factor matrices except $\mathbf{U}_{n}$. The MTTKRP is followed by the Hadamard product of the Gramians of each factor matrix (${\mathbf{U}_{i}}^{T}\mathbf{U}_{i}$) in line 1. Factor matrix $\mathbf{U}_{n}$ is updated by solving the linear system in line 1. At the completion of an iteration, i.e., a full pass over all $N$ modes, the error between the model and the target tensor is computed (line 1) using the efficient formula derived in (Phan et al., 2013).

Assuming a small number of components ($R$), the most expensive step is the MTTKRP (line 1). This step involves $2R\prod_{i}I_{i}$ FLOPs (ignoring, for the sake of simplicity, the lower order of FLOPs required for the computation of the KRP). The operation touches slightly more than $\prod_{i}I_{i}$ memory locations, resulting in an arithmetic intensity less than $2R$ FLOPs per memory reference. Thus, unless $R$ is sufficiently large, the speed of the computation will be limited by the memory bandwidth rather than the speed of the processor. The CP-ALS algorithm is inherently memory-bound for small $R$, regardless of how it is implemented.

The impact on performance of the memory-bound nature of MTTKRP is demonstrated in Fig. 1, which shows the computational efficiency of a particular implementation of MTTKRP as a function of the number of components (for a tensor of size $50\times 200\times 200$). Efficiency is defined as the ratio of the performance achieved by MTTKRP (in FLOPs/sec), relative to the Theoretical Peak Performance (TPP, see below) of the machine, i.e.,

The TPP of a machine is defined as the maximum number of (double precision) floating point operations the machine can perform in one second. Table 1 shows the TPP for our particular machine (see Sec. 5 for details).

In Fig. 1, we see that the efficiency of MTTKRP tends to increase with the number of components, $R$, until eventually reaching a plateau. On this machine, the plateau is $R\geq 60$ at $\approx 70\%$ efficiency for one thread and $R\geq 300$ at $\approx 35\%$ efficiency for 24 threads. For $R\leq 20$, which is common in applications, the efficiency is well below the TPP.

When fitting multiple CP models to the same underlying tensor, the Concurrent ALS (CALS) technique can improve the efficiency if the number of components is not large enough for CP-ALS to reach its performance plateau (Psarras et al., 2020). A need to fit multiple models to the same tensor arises, for example, when trying different initial guesses or when trying different numbers of components.

The gist of CALS can be summarized as follows (see (Psarras et al., 2020) for details). Suppose $K$ independent instances of CP-ALS have to be executed on the same underlying tensor. Rather than running them sequentially or in parallel, run them in lock-step fashion as follows. Advance every CP-ALS process one iteration before proceeding to the next iteration. One CALS iteration entails $K$ CP-ALS iterations (one iteration per model). Each CP-ALS iteration in turn contains one MTTKRP operation, so one CALS iteration also entails $K$ MTTKRP operations. But these MTTKRPs all involve the same tensor and can therefore be fused into one bigger MTTKRP operation (see Eq. 3 of (Psarras et al., 2020)). The performance of the fused MTTKRP depends on the sum total of components, i.e., $\sum_{i=1}^{K}R_{i}$, where $R_{i}$ is the number of components in model $i$. Due to the performance profile of MTTKRP (see Fig. 1), the fused MTTKRP is expected to be more efficient than each of the $K$ smaller operations it replaces.

The following example illustrates the impact on efficiency of MTTKRP fusion. Given a target tensor of size $50\times 200\times 200$, $K=50$ models to fit, and $R_{i}=5$ components in each model, the efficiency for each of the $K$ MTTKRPs in CP-ALS is about $15\%$ ($3\%$) for 1 (24) threads (see Fig.1). The efficiency of the fused MTTKRP in CALS will be as observed for $R=\sum_{i=1}^{K}R_{i}=250$, i.e., $60\%$ ($30\%$) for 1 (24) threads. Since the MTTKRP operation dominates the cost, CALS is expected to be $\approx 4\times$ ($\approx 10\times$) faster than CP-ALS for 1 (24) threads.

Algorithm 2 shows a baseline (inefficient) application of leave-one-out jackknife resampling to a CP model. For details, see (Riu and Bro, 2003). The inputs are a target tensor $\boldsymbol{\mathcal{T}}$, an overall CP model $P$ fitted to all of $\boldsymbol{\mathcal{T}}$, and a sampled mode $\hat{n}\in\{1,2,\ldots,N\}$. For each sample $p\in\{1,2,\ldots,I_{\hat{n}}\}$, the algorihm removes the slice corresponding to the sample from tensor $\boldsymbol{\mathcal{T}}$ (line 2) and model $P$ (line 2) and fits a reduced model $P_{-p}$ (lines 2–2) to the reduced tensor $\boldsymbol{\mathcal{\hat{T}}}$ using regular CP-ALS. After fitting all submodels, the standard deviation of every factor matrix except $\mathbf{U}_{\hat{n}}$ is computed from the $I_{\hat{n}}$ submodels in $\mathcal{Q}$ (line 2). The only costly part of Algorithm 2 is the repeated calls to CP-ALS in line 2.

Let $\boldsymbol{\mathcal{T}}$ be an $N$-mode tensor with a corresponding CP model $\mathbf{A}_{1},\ldots,\mathbf{A}_{N}$. Let $\boldsymbol{\mathcal{\hat{T}}}$ be identical to $\boldsymbol{\mathcal{T}}$ except for one sample (with index $p$) removed from the sampled mode $\hat{n}\in\{1,2,\ldots,N\}$. Let $\mathbf{\hat{A}}_{1},\ldots,\mathbf{\hat{A}}_{N}$ be the CP submodel corresponding to the resampled tensor $\boldsymbol{\mathcal{T}}$.

When fitting a CP model to $\boldsymbol{\mathcal{T}}$ using CP-ALS, the MTTKRP for mode $n$ is given by

Similarly, when fitting a model to $\boldsymbol{\mathcal{\hat{T}}}$, the MTTKRP for mode $n$ is given by

Can $\mathbf{\hat{M}}_{n}$ be computed from $\mathbf{T}_{(n)}$ instead of $\mathbf{\hat{T}}_{(n)}$? As we will see, the answer is yes. We separate two cases: $n=\hat{n}$ and $n\neq\hat{n}$.

Case I: $n=\hat{n}$. The slice of $\boldsymbol{\mathcal{T}}$ removed when resampling corresponds to a row of the unfolding $\mathbf{T}_{(n)}=\mathbf{T}_{(\hat{n})}$. To see this, note that element $\boldsymbol{\mathcal{T}}(i_{1},i_{2},\dots,i_{N})$ corresponds to element $\mathbf{T}_{(n)}(i_{n},j)$ of its mode-$n$ unfolding (Kolda and Bader, 2009), where

When we remove sample $p$, then $\mathbf{\hat{T}}_{(n)}$ will be identical to $\mathbf{T}_{(n)}$ except that row $p$ from the latter is missing in the former. In other words, $\mathbf{\hat{T}}_{(n)}=\mathbf{E}_{p}\mathbf{T}_{(n)}$, where $\mathbf{E}_{p}$ is the matrix that removes row $p$. We can therefore compute $\mathbf{\hat{M}}_{n}$ by replacing $\mathbf{\hat{T}}_{(n)}$ with $\mathbf{T}_{(n)}$ in (2) and then discarding row $p$ from the result:

Case II: $n\neq\hat{n}$. The slice of $\boldsymbol{\mathcal{T}}$ removed when resampling corresponds to a set of columns in the unfolding $\mathbf{T}_{(n)}$. One could in principle remove these columns to obtain $\mathbf{\hat{T}}_{(n)}$. But instead of explicitly removing sample $p$ from $\boldsymbol{\mathcal{T}}$, we can simply zero out the corresponding slice of $\boldsymbol{\mathcal{T}}$. To give the CP model matching dimensions, we need only insert a row of zeros at index $p$ in factor matrix $\hat{n}$. Crucially, the zeroing out of slice $p$ is superfluous. In the MTTKRP, the elements that should have been zeroed out will be multiplied with zeros in the Khatri-Rao product generated by the row of zeros insert in factor matrix $\hat{n}$. Thus, to compute $\mathbf{\hat{M}}_{n}$ in (2) we (a) replace $\mathbf{\hat{T}}_{(n)}$ with $\mathbf{T}_{(n)}$ and (b) insert a row of zeros at index $p$ in factor matrix $\mathbf{\hat{A}}_{\hat{n}}$.

In summary, we have shown that it is possible to compute $\mathbf{\hat{M}}_{n}$ in (2) without referencing the reduced tensor. There is an overhead associate with extra arithmetic. For the case $n=\hat{n}$, we compute numbers that are later discarded. For the case $n\neq\hat{n}$, we do some arithmetic with zeros.

Based on the observations above, the CALS algorithm (Psarras et al., 2020) can be modified to facilitate the concurrent fitting of all jackknife submodels. Algorithm 3 incorporates the necessary changes. In the end, extending CALS to support jackknife comes down to these localized changes (colored red in Algorithm 3):

• •

  Insert a row of zeros in one of the factor matrices,

• •

  periodically zero out the padded row to keep it zero, and

• •

  adjust the error formula to compute the submodel error.

We remark that JK-CALS can be straightforwardly extended to delete-$d$ jackknife. Instead of padding and periodically zeroing out one row, we pad and periodically zero out $d$ rows.

While Algorithm 3 benefits from improved MTTKRP efficiency, the padding results in extra arithmetic operations. Let $d$ denote the number of removed samples ($d=1$ corresponds to leave-one-out). For the sake of simplicity, assume that the integer $d$ divides $I_{\hat{n}}$. There are $I_{\hat{n}}/d$ submodels, each with $R$ components. The only costly part is the MTTKRP.

The MTTKRPs in JK-ALS (for all submodels combined) requires

FLOPs. Meanwhile, the fused MTTKRP in JK-CALS requires

FLOPs. The ratio of the latter to the former comes down to

since $d\leq I_{\hat{n}}/2$ in delete-$d$ jackknife. Thus, in the worst case, JK-CALS requires less than twice the FLOPs of JK-ALS. More typically, the overhead is negligible.

In this first experiment, we use three synthetic tensors of size $50\times m\times m$ with $m\in\{100,200,400\}$, referred to as “small”, “medium” and “large” tensors, respectively. The samples are in the first mode. The other modes contain variables. The number of samples is kept low, since leave-one-out jackknife is usually performed on a small number of samples (usually $<100$), while there can be arbitrarily many variables.

For each tensor, we perform jackknife on four models with varying number of components ($R\in\{3,5,7,9\}$). This range of component counts is typical in applications. In practice, it is often the case that multiple models are fitted to the target tensor, and many of those models are then further analyzed using jackknife. For this reason, we perform jackknife on each model individually, as well as to all models simultaneously (denoted by “All” in the figures), to better simulate multiple real-world application scenarios. In this experiment, the termination criteria based on maximum number of iterations and tolerance are ignored; instead, all models are forced to go through exactly 100 iterations, typically a small number of iterations for small values of tolerance (i.e., most models require more than 100 iterations). The reason for this choice is that we aim to isolate the performance difference of the methods tested; therefore, we maintain a consistent amount of workload throughout the experiment. (Tolerance and maximum number of iterations are instead used later on in the application experiments.)

For comparison, we perform jackknife using three methods: JK-ALS, JK-OALS and JK-CALS. JK-OALS uses OpenMP to take advantage of the inherent parallelism when fitting multiple submodels. This method is only used for multi-threaded and GPU experiments, and we are only going to focus on its performance, ignoring the memory overhead associated with it.

Fig. 2 shows results for single threaded execution; in this case, JK-OALS is absent. JK-CALS consistently outperforms JK-ALS for all tensor sizes and workloads. Specifically, for any fixed amount of workload—i.e., a model of a specific number of components—JK-CALS exhibits increasing speedups compared to JK-ALS, as the tensor size increases. For example, for a model with 5 components, JK-CALS is $2.9$, $3$, $5.2$ times faster than JK-ALS, for the small, medium and large tensor sizes, respectively.

Fig. 3 shows results for multi-threaded execution, using 24 threads. In this case, JK-CALS outperforms the other two implementations (JK-ALS and JK-OALS) for the medium and large tensors, for all workloads (number of components), exhibiting speedups up to $35\times$ and $8\times$ compared to JK-ALS and JK-OALS, respectively. For the small tensor ($50\times 100\times 100$) and small workloads ($R\leq 7$), JK-CALS is outperformed by JK-OALS; for $R=3$, it is also outperformed by JK-ALS. Investigating this further, for the small tensor and $R=3$ and $5$, the parallel speedup (the ratio between single threaded and multi-threaded execution time) of JK-CALS is $0.3\times$ and $0.7\times$ for 24 threads. However, for $12$ threads, the corresponding timings are $0.28$ and $0.27$ seconds, resulting in speedups of $2.7\times$ and $3.7\times$ respectively. This points to two main reasons for the observed performance of JK-CALS in these cases: a) the amount of available computational resources (24 threads) is disproportionately high compared to the volume of computation to be performed and b) because of the small amount of overall computation, the small overhead associated with the CALS methodology becomes more significant.

That being said, even for the small tensor, as the amount of workload increases—already for a single model with 9 components—JK-CALS again becomes the fastest method. Finally, similarly to the single threaded case, as the size of the tensor increases, so do the speedups achieved by JK-CALS over the other two methods.

Fig. 4 shows results when the GPU is used to perform MTTKRP for all three methods; in this case, all 24 threads are used on the CPU. For the small tensor and small workloads ($R\leq 5$), there is not enough computation to warrant the shipping of data to and from the GPU, resulting in higher execution times compared to multi-threaded execution; for all other cases, all methods have reduced execution time when using the GPU compared to the execution on 24 threads. Furthermore, in those cases, JK-CALS is consistently faster than its counterparts, exhibiting the largest speedups when the workload is at its highest (“All”), with values of $10\times$, $7\times$, $7\times$ compared to JK-OALS, and $12\times$, $10\times$, $9\times$ compared to JK-ALS, for the small, medium and large tensors, respectively.

In this second experiment, we selected two tensors of size $89\times 97\times 549$ and $44\times 2700\times 200$ from the field of Chemometrics (Acar et al., 2008; Skov et al., 2008). In this field it is common to fit multiple, randomly initialized models in a range of low components (e.g. $R\in\{1,2,\ldots,20\}$, $10$–$20$ models for each $R$, and then analyze (e.g., using jackknife) those models that might be of particular interest (often those with components close to the expected rank of the target tensor); in the tensors we consider, the expected rank is $5$ and $20$, respectively. To mimic the typical workflow of practitioners, we fitted three models to each tensor, of components $R\in\{4,5,6\}$ and $R\in\{19,20,21\}$, respectively, and used the three methods (JK-ALS, JK-OALS and JK-CALS) to apply jackknife resampling to the fitted models. The values for tolerance and maximum number of iterations were set according to typical values for the particular field, namely $10^{-6}$ and $1000$, respectively.

In Fig. 5 we report the execution time for 1 thread, 24 threads, and GPU + 24 threads. For both datasets and for all configurations, JK-CALS is faster than the other two methods. Specifically, when compared to JK-ALS over the two tensors, JK-CALS achieves speedups of $5.4\times$ and $2\times$ for single threaded execution, $10\times$ and $2.8\times$ for 24-threaded execution. Similarly, when compared to JK-OALS, JK-CALS achieves speedups of $2.7\times$ and $4.8\times$ for 24-threaded execution. Finally, JK-CALS takes advantage of the GPU the most, exhibiting speedups of $17.5\times$ and $3.7\times$ over JK-ALS, and $9\times$ and $2\times$ over JK-OALS, for GPU execution.