TuckerMPI: A Parallel C++/MPI Software Package for Large-scale Data Compression via the Tucker Tensor Decomposition

TL;DR

TuckerMPI is a parallel C++/MPI software that efficiently compresses massive tensor-structured data from high-fidelity simulations using Tucker decomposition, achieving high compression ratios and enabling partial data reconstruction on single nodes.

Contribution

This work introduces TuckerMPI, a novel parallel software for tensor-based data compression that handles multi-terabyte datasets with high efficiency and supports partial data reconstruction.

Findings

Achieved 100-200,000× compression ratios on multi-terabyte datasets.

Compressed data in less time than reading from parallel filesystems.

Enabled partial data reconstruction on single nodes for visualization and analysis.

Abstract

Our goal is compression of massive-scale grid-structured data, such as the multi-terabyte output of a high-fidelity computational simulation. For such data sets, we have developed a new software package called TuckerMPI, a parallel C++/MPI software package for compressing distributed data. The approach is based on treating the data as a tensor, i.e., a multidimensional array, and computing its truncated Tucker decomposition, a higher-order analogue to the truncated singular value decomposition of a matrix. The result is a low-rank approximation of the original tensor-structured data. Compression efficiency is achieved by detecting latent global structure within the data, which we contrast to most compression methods that are focused on local structure. In this work, we describe TuckerMPI, our implementation of the truncated Tucker decomposition, including details of the data distribution and in-memory layouts, the parallel and serial implementations of the key kernels, and analysis of the storage, communication, and computational costs. We test the software on 4.5 terabyte and 6.7 terabyte data sets distributed across 100s of nodes (1000s of MPI processes), achieving compression rates between 100-200,000$\times$ which equates to 99-99.999% compression (depending on the desired accuracy) in substantially less time than it would take to even read the same dataset from a parallel filesystem. Moreover, we show that our method also allows for reconstruction of partial or down-sampled data on a single node, without a parallel computer so long as the reconstructed portion is small enough to fit on a single machine, e.g., in the instance of reconstructing/visualizing a single down-sampled time step or computing summary statistics.

Tables

Table 1. Table 1. Local tensor sizes for a 3 × 4 × 3 × 2 3 4 3 2 3\times 4\times 3\times 2 tensor distributed on a 2 × 2 × 2 × 1 2 2 2 1 2\times 2\times 2\times 1 processor grid, demonstrating how uneven sizes are handled.

Proc. ID	Proc. Linear ID	Local Size
$(p_0,p_1,p_2,p_3)$	$p^{\prime }$	${\overline{J}}_0\times {\overline{J}}_1\times {\overline{J}}_2\times {\overline{J}}_3$
$(0,0,0,0)$	0	$2\times 2\times 2\times 2$
$(1,0,0,0)$	1	$1\times 2\times 2\times 2$
$(0,1,0,0)$	2	$2\times 2\times 2\times 2$
$(1,1,0,0)$	3	$1\times 2\times 2\times 2$
$(0,0,1,0)$	4	$2\times 2\times 1\times 2$
$(1,0,1,0)$	5	$1\times 2\times 1\times 2$
$(0,1,1,0)$	6	$2\times 2\times 1\times 2$
$(1,1,1,0)$	7	$1\times 2\times 1\times 2$

Table 2. Table 2. Communication costs in MPI model, where W 𝑊 W is the local input data size and P 𝑃 P is the number of processors.

Send/Receive	$\alpha +\beta W$
Reduce/All-Reduce	$2\alpha \log P+\left(2\beta +\gamma \right)\frac{P-1}{P}W$
Reduce-Scatter	$\alpha \log P+\left(\beta +\gamma \right)\frac{P-1}{P}W$
All-to-All	$\alpha \left(P-1\right)+\beta \frac{P-1}{P}W$

Table 3. Table 3. Asymptotic costs of algorithms for TTM and Gram with respect to mode n 𝑛 n for a tensor with global dimensions J 0 × J 1 × ⋯ × J N − 1 superscript subscript 𝐽 0 absent superscript subscript 𝐽 1 absent ⋯ superscript subscript 𝐽 𝑁 1 absent {J}_{0}^{\phantom{{}^{\prime}}}\times{J}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{J}_{N-1}^{\phantom{{}^{\prime}}} and processor grid with dimensions P 0 × P 1 × ⋯ × P N − 1 superscript subscript 𝑃 0 absent superscript subscript 𝑃 1 absent ⋯ superscript subscript 𝑃 𝑁 1 absent {P}_{0}^{\phantom{{}^{\prime}}}\times{P}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{P}_{N-1}^{\phantom{{}^{\prime}}} . We omit leading constant factors and lower order terms. The columns correspond to per-processor costs: number of floating point operations, number of words communicated, number of messages communicated, and amount of temporary local memory required, respectively. Recall J \oast = ∏ n ∈ [ N ] J n superscript 𝐽 \oast subscript product 𝑛 delimited-[] 𝑁 subscript 𝐽 𝑛 J^{\oast}=\prod_{n\in[\,N\,]}J_{n} and J n ⊘ = J \oast / J n superscript subscript 𝐽 𝑛 ⊘ superscript 𝐽 \oast subscript 𝐽 𝑛 {J}_{n}^{\oslash}=J^{\oast}/J_{n} , with analogous definitions for P \oast superscript 𝑃 \oast P^{\oast} and P n ⊘ superscript subscript 𝑃 𝑛 ⊘ {P}_{n}^{\oslash} .

	Flops	Bandwidth	Latency	Memory
TTM (Reduce-Scatter)	$\frac{J^{\text{oast}}}{P^{\text{oast}}}K$	$\frac{J_n^{\oslash }\mathrm{missing}}{P_n^{\oslash }\mathrm{missing}}K$	$\log P_n$	$\frac{J_n^{\oslash }\mathrm{missing}}{P_n^{\oslash }\mathrm{missing}}K$
TTM (Reduce)	$\frac{J^{\text{oast}}}{P^{\text{oast}}}K$	$\frac{J_n^{\oslash }\mathrm{missing}}{P_n^{\oslash }\mathrm{missing}}K$	$P_n\log P_n$	$\frac{J_n^{\oslash }\mathrm{missing}}{P^{\text{oast}}}K$
Gram (Original)	$\frac{J^{\text{oast}}}{P^{\text{oast}}}{J}_n$	$\frac{J^{\text{oast}}}{P^{\text{oast}}}{P}_n+\frac{J_n^2}{P_n}$	$P_n$	$\frac{J^{\text{oast}}}{P^{\text{oast}}}+{\left(\frac{J_n}{P_n}\right)}^2$
Gram (New)	$\frac{J^{\text{oast}}}{P^{\text{oast}}}{J}_n$	$\frac{J^{\text{oast}}}{P^{\text{oast}}}+J_n^2$	$P_n$	$\frac{J^{\text{oast}}}{P^{\text{oast}}}+J_n^2$

Table 4. Table 4. Experimental setup

Dataset	Overall	Total	Number of		Storage per
Name	Tensor Size	Storage	Nodes	Processes	Node	Process
SP	$500\times 500\times 500\times 11\times 400$	4.4 TB	250	4000	17.6 GB	1.1 GB
JICF	$1500\times 2080\times 1500\times 18\times 10$	6.7 TB	350	5600	19.3 GB	1.2 GB

Table 5. (a) Data sets to be compressed and the number of parallel nodes to be used in the compression experiments.

Dataset	Overall	Total	Number of		Storage per
Name	Tensor Size	Storage	Nodes	Processes	Node	Process
SP	$500\times 500\times 500\times 11\times 400$	4.4 TB	250	4000	17.6 GB	1.1 GB
JICF	$1500\times 2080\times 1500\times 18\times 10$	6.7 TB	350	5600	19.3 GB	1.2 GB

Table 6. (b) Three different processor grid configurations per dataset, to test the efficiency of the compression. The local tensor size may vary, but here we list the largest local size in each dimension.

Dataset	Proc. Config. Name	Processor Grid	Local Tensor Size
SP	A	$1\times 1\times 40\times 1\times 100$	$500\times 500\times 13\times 11\times 4$
	B	$10\times 8\times 5\times 1\times 10$	$50\times 63\times 100\times 11\times 40$
	C	$40\times 10\times 1\times 1\times 10$	$13\times 50\times 500\times 11\times 40$
JICF	A	$1\times 16\times 35\times 1\times 10$	$1500\times 130\times 43\times 18\times 1$
	B	$10\times 8\times 7\times 1\times 10$	$150\times 260\times 215\times 18\times 1$
	C	$35\times 16\times 1\times 1\times 10$	$43\times 130\times 1500\times 18\times 1$

Table 7. Table 5. Intermediate tensor sizes in the partial reconstruction for two different TTM orderings, using the high compression scenario ( ϵ italic-ϵ \epsilon =1e-2) for the SP dataset. The TTM ordering can make a dramatic different in memory usage.

	Ordering 01234		Ordering 43120
TTM	Dimensions	Intermed. Size (GB)	Dimensions	Intermed. Size (GB)
None	$30\times 38\times 35\times 6\times 11$	0.02	$30\times 38\times 35\times 6\times 11$	0.02
1st	$500\times 38\times 35\times 6\times 11$	0.35	$30\times 38\times 35\times 6\times 1$	0.01
2nd	$500\times 500\times 35\times 6\times 11$	4.62	$30\times 38\times 35\times 1\times 1$	0.01
3rd	$500\times 500\times 500\times 6\times 11$	66.0	$30\times 500\times 35\times 1\times 1$	0.01
4th	$500\times 500\times 500\times 1\times 11$	11.0	$30\times 500\times 500\times 1\times 1$	0.06
5th	$500\times 500\times 500\times 1\times 1$	1.00	$500\times 500\times 500\times 1\times 1$	1.00

Table 8. Table 6. Partial reconstruction results on the SP dataset with a single MPI process, computing a partial reconstruction of size 500 × 500 × 500 × 1 × 1 500 500 500 1 1 500\times 500\times 500\times 1\times 1 (1 GB).

Compression	Mode	Max Memory	Compute	I/O Time (sec)
Scenario	Order	Usage (GB)	Time (sec)	Input	Output
High $ϵ$=1e-2	01234	out of memory		0.05 (21.5 MB)	1.01 (1 GB)
	03421	1.08	1.55
	24013	7.00	8.88
	34021	7.00	1.19
	43120	1.08	1.05
Low $ϵ$=1e-4	01234	out of memory		12.23 (10.7 GB)
	03421	out of memory
	24013	53.62	129.70
	34021	12.26	5.83
	43120	10.81	4.72

Table 9. Table 7. Summary statistic results on the SP dataset with a single MPI process, computing the average fraction of carbon dioxide across all time steps and entire physical grid. Relative error is in comparison to the value computed based on the original data.

Compression	Mode	Max Memory	Compute	I/O Time (sec)		Relative
Scenario	Order	Usage (GB)	Time (sec)	Input	Output	Error
High ($ϵ$=1e-2)	01234	0.022	0.16	0.69 (21.5 MB)	0.04 (8 B)	6.7e-6
Low ($ϵ$=1e-4)	01234	10.84	2.63	13.02 (10.7 GB)		5.7e-8

Table 10. Table 8. Full reconstruction results of the SP dataset using 250 nodes/4000 MPI processes. The processor grid is size 1 × 1 × 40 × 1 × 100 1 1 40 1 100 1\times 1\times 40\times 1\times 100 . The “Max Memory Usage” is per process.

Compression	Mode	Max Memory	Compute	I/O Time (sec)
Scenario	Order	Usage (GB)	Time (sec)	Input	Output
High $ϵ$=1e-2	12034	2.57	320.34	4.02 (21.5 MB)	3118.57 (4.4 TB)
	34120	1.21	27.51
	41203	1.77	16.39
	42103	1.77	5.04
Low $ϵ$=1e-4	12034	out of memory		15.24 (10.7 GB)
	34120	1.37	82.46
	41203	1.88	53.74
	42103	1.88	19.43

Equations

$$\Vert \mathcal{X}-\hat{\mathcal{X}}\Vert \le ϵ\Vert \mathcal{X}\Vert .$$

$$I^{\text{oast}}=\prod _{n\in [N]}{I}_n.$$

$$I_n^{\oslash }$$

$$I_n^{\text{varolessthan}}$$

$$I_n^{\text{varogreaterthan}}$$

$$i_n^{\prime }=\sum _{k<n}{i}_k\cdot I_k^{\text{varolessthan}}+\sum _{k>n}{i}_k\cdot I_k^{\text{varolessthan}}/I_n.$$

$$\mathcal{Y}=\mathcal{X}{\times }_n\mathbf{U}\iff {\mathbf{Y}}_{(n)}=\mathbf{U}{\mathbf{X}}_{(n)}.$$

$$\mathcal{X}{\times }_m\mathbf{U}{\times }_n\mathbf{V}=\mathcal{X}{\times }_n\mathbf{V}{\times }_m\mathbf{U}\text{ for }m\ne n.$$

$$\mathcal{X}{\times }_n\mathbf{U}{\times }_n\mathbf{V}=\mathcal{X}{\times }_n\mathbf{VU},$$

$$\mathcal{X}\approx \hat{\mathcal{X}}\equiv \mathcal{G}{\times }_0{\mathbf{U}}_0{\times }_1{\mathbf{U}}_1\cdots {\times }_{N-1}{\mathbf{U}}_{N-1}.$$

$$\mathcal{G}=\mathcal{X}{\times }_0{\mathbf{U}}_0^{\mathsf{T}}{\times }_1{\mathbf{U}}_1^{\mathsf{T}}\cdots {\times }_{N-1}{\mathbf{U}}_{N-1}^{\mathsf{T}}.$$

$$\hat{\mathcal{X}}=\mathcal{X}{\times }_0{\mathbf{U}}_0^{}{\mathbf{U}}_0^{\mathsf{T}}{\times }_1{\mathbf{U}}_1^{}{\mathbf{U}}_1^{\mathsf{T}}\cdots {\times }_{N-1}{\mathbf{U}}_{N-1}^{}{\mathbf{U}}_{N-1}^{\mathsf{T}}.$$

$${\tilde{\mathcal{G}}}_n=\mathcal{X}{\times }_0{\mathbf{U}}_0^{\mathsf{T}}\cdots {\times }_n{\mathbf{U}}_n^{\mathsf{T}}\text{for}n\in [N].$$

$${\tilde{\mathcal{X}}}_n=\mathcal{G}{\times }_0{\mathbf{U}}_0\cdots {\times }_n{\mathbf{U}}_n=\mathcal{X}{\times }_0{\mathbf{U}}_0^{}{\mathbf{U}}_0^{\mathsf{T}}\cdots {\times }_n{\mathbf{U}}_n^{}{\mathbf{U}}_n^{\mathsf{T}}\text{for}n\in [N].$$

$$\Bigl{\|}\;{\bm{{\mathscr{{X}}}}}-{\bm{\hat{\mathscr{{X}}}}}\;\Bigr{\|}^{2}=\Bigl{\|}\;{\bm{{\mathscr{{X}}}}}\times_{0}\Bigl{(}{\bm{\mathbf{{I}}}}-{\mathbf{U}}_{0}^{\phantom{T}}\!{\mathbf{U}}_{0}^{\sf T}\Bigr{)}\;\Bigr{\|}^{2}+\Bigl{\|}\;{\bm{\tilde{\mathscr{{X}}}}}_{0}\times_{1}\Bigl{(}{\bm{\mathbf{{I}}}}-{\mathbf{U}}_{1}^{\phantom{T}}\!{\mathbf{U}}_{1}^{\sf T}\Bigr{)}\;\Bigr{\|}^{2}+\cdots+\Bigl{\|}\;{\bm{\tilde{\mathscr{{X}}}}}_{N-2}\times_{N-1}\Bigl{(}{\bm{\mathbf{{I}}}}-{\mathbf{U}}_{N-1}^{\phantom{T}}\!{\mathbf{U}}_{N-1}^{\sf T}\Bigr{)}\;\Bigr{\|}^{2},$$

$$\Bigl{\|}\;{\bm{{\mathscr{{X}}}}}-{\bm{\hat{\mathscr{{X}}}}}\;\Bigr{\|}^{2}\leq\sum_{n\in[\,N\,]}\Bigl{\|}\;{\bm{{\mathscr{{X}}}}}\times_{n}\Bigl{(}{\bm{\mathbf{{I}}}}-{\mathbf{U}}_{n}^{\phantom{T}}\!{\mathbf{U}}_{n}^{\sf T}\Bigr{)}\;\Bigr{\|}^{2}.$$

$$\frac{\Vert \mathcal{X}-\hat{\mathcal{X}}\Vert }{\Vert \mathcal{X}\Vert }\le ϵ.$$

$$\mathbf{S}\equiv {\mathbf{X}}_{(n)}^{}{\mathbf{X}}_{(n)}^{\mathsf{T}}=\mathbf{V}\mathbf{\Lambda }{\mathbf{V}}^{\mathsf{T}}.$$

$${\mathbf{U}}_n=\mathbf{V}(:,1:R_n)\text{where}{R}_n=\underset{R\in [I_n]}{\min }R\text{subject to }\sum _{i=R+1}^{I_n}{\lambda }_i\le {ϵ}^2\Vert \mathcal{X}{\Vert }^2/N.$$

$$\Bigl{\|}\;{\bm{{\mathscr{{X}}}}}\times_{n}\Bigl{(}{\bm{\mathbf{{I}}}}-{\mathbf{U}}_{n}^{\phantom{T}}\!{\mathbf{U}}_{n}^{\sf T}\Bigr{)}\;\Bigr{\|}^{2}\leq\epsilon^{2}\|{\bm{{\mathscr{{X}}}}}\|^{2}/N.$$

$$\Bigl{\|}\;{\bm{\tilde{\mathscr{{X}}}}}_{n-1}\times_{n}\Bigl{(}{\bm{\mathbf{{I}}}}-{\mathbf{U}}_{n}^{\phantom{T}}\!{\mathbf{U}}_{n}^{\sf T}\Bigr{)}\;\Bigr{\|}^{2}=\Bigl{\|}\;{\bm{\tilde{\mathscr{{G}}}}}_{n-1}\times_{n}\Bigl{(}{\bm{\mathbf{{I}}}}-{\mathbf{U}}_{n}^{\phantom{T}}\!{\mathbf{U}}_{n}^{\sf T}\Bigr{)}\;\Bigr{\|}^{2}.$$

$$\Bigl{\|}\;{\bm{{\mathscr{{X}}}}}-{\bm{\hat{\mathscr{{X}}}}}\;\Bigr{\|}^{2}=\Bigl{\|}\;{\bm{{\mathscr{{X}}}}}\times_{0}\Bigl{(}{\bm{\mathbf{{I}}}}-{\mathbf{U}}_{0}^{\phantom{T}}\!{\mathbf{U}}_{0}^{\sf T}\Bigr{)}\;\Bigr{\|}^{2}+\Bigl{\|}\;{\bm{\tilde{\mathscr{{G}}}}}_{0}\times_{1}\Bigl{(}{\bm{\mathbf{{I}}}}-{\mathbf{U}}_{1}^{\phantom{T}}\!{\mathbf{U}}_{1}^{\sf T}\Bigr{)}\;\Bigr{\|}^{2}+\cdots+\Bigl{\|}\;{\bm{\tilde{\mathscr{{G}}}}}_{N-2}\times_{N-1}\Bigl{(}{\bm{\mathbf{{I}}}}-{\mathbf{U}}_{N-1}^{\phantom{T}}\!{\mathbf{U}}_{N-1}^{\sf T}\Bigr{)}\;\Bigr{\|}^{2},$$

$$\mathbf{S}\equiv {\mathbf{Y}}_{(n)}^{}{\mathbf{Y}}_{(n)}^{\mathsf{T}}=\mathbf{V}\mathbf{\Lambda }{\mathbf{V}}^{\mathsf{T}}\text{where}\mathcal{Y}={\tilde{\mathcal{G}}}_{n-1}.$$

$$\Bigl{\|}\;{\bm{{\mathscr{{X}}}}}-{\bm{\hat{\mathscr{{X}}}}}_{\rm{HOSVD}}\;\Bigr{\|}\leq\sqrt{N}\;\Bigl{\|}\;{\bm{{\mathscr{{X}}}}}-{\bm{\hat{\mathscr{{X}}}}}_{\rm{opt}}\;\Bigr{\|}$$

$$\Bigl{\|}\;{\bm{{\mathscr{{X}}}}}-{\bm{\hat{\mathscr{{X}}}}}_{\rm{ST-HOSVD}}\;\Bigr{\|}\leq\sqrt{N}\;\Bigl{\|}\;{\bm{{\mathscr{{X}}}}}-{\bm{\hat{\mathscr{{X}}}}}_{\rm{opt}}\;\Bigr{\|}$$

$$j^{\prime}=\operatorname{idx2lin}\bigl{\{}\,(j_{0},j_{1},\dots,j_{N-1})\,,\,(J_{0},J_{1},\dots,J_{N-1})\,\bigr{\}}\equiv\sum_{n\in[\,N\,]}j_{n}\cdot{J}_{n}^{\varolessthan},$$

$$(j_{0},j_{1},\dots,j_{N-1})=\operatorname{lin2idx}\bigl{\{}\,j^{\prime}\,,\,(J_{0},J_{1},\dots,J_{N-1})\,\bigr{\}}.$$

$${\bm{\mathbf{{Y}}}}_{(0)}=\left[\begin{array}[]{*{23}{r|}r}0&3&6&9&12&15&18&21&24&27&30&33&36&39&42&45&48&51&54&57&60&63&66&69\\ 1&4&7&10&13&16&19&22&25&28&31&34&37&40&43&46&49&52&55&58&61&64&67&70\\ 2&5&8&11&14&17&20&23&26&29&32&35&38&41&44&47&50&53&56&59&62&65&68&71\\ \end{array}\right].$$

$${\bm{\mathbf{{Y}}}}_{(1)}=\left[\begin{array}[]{*{5}{rrr|}rrr}0&1&2&12&13&14&24&25&26&36&37&38&48&49&50&60&61&62\\ 3&4&5&15&16&17&27&28&29&39&40&41&51&52&53&63&64&65\\ 6&7&8&18&19&20&30&31&32&42&43&44&54&55&56&66&67&68\\ 9&10&11&21&22&23&33&34&35&45&46&47&57&58&59&69&70&71\\ \end{array}\right].$$

$${\bm{\mathbf{{Y}}}}_{(2)}=\left[\begin{array}[]{rrrrrrrrrrrr|rrrrrrrrrrrr}0&1&2&3&4&5&6&7&8&9&10&11&36&37&38&39&40&41&42&43&44&45&46&47\\ 12&13&14&15&16&17&18&19&20&21&22&23&48&49&50&51&52&53&54&55&56&57&58&59\\ 24&25&26&27&28&29&30&31&32&33&34&35&60&61&62&63&64&65&66&67&68&69&70&71\\ \end{array}\right].$$

$${\bm{\mathbf{{Y}}}}_{(3)}=\left[\begin{array}[]{rrrrrrrrrrrrrrrrrrrrrrr}0&1&2&3&4&5&6&7&8&9&10&11&\cdots&26&27&28&29&30&31&32&33&34&35\\ 36&37&38&39&40&41&42&43&44&45&46&47&\cdots&62&63&64&65&66&67&68&69&70&71\\ \end{array}\right].$$

Full text

Data compression, Distributed algorithms

Tucker decomposition, tensor decomposition, higher-order singular value decomposition (HOSVD)

††copyright: none††journal: TOMS

1. Introduction

The convergence of ever-faster computational platforms and algorithmic advancements in scientific simulations have led to a data deluge — it is now possible to produce very fine-grained and detailed simulations, resulting in terabyte-sized or larger datasets. Consider the problem we discuss later on in this work: a simulation on a three-dimensional rectangular grid of size $500\times 500\times 500$, tracking 11 variables for 400 time steps. Even this modest-sized simulation yields 4 TB of data in double precision. With current computational platforms, this data cannot easily be stored, moved, visualized, or analyzed. Nevertheless, it is well known to simulation scientists that their massive datasets have extensive latent structure and are therefore highly compressible. The problem is how to discover the redundancies automatically.

Most compression methods focus on compressing local structure with very little loss in precision. Fout, Ma, and Ahrens (FoMaAh05, ) do multivariate volume block data reduction to take advantage of local multiway structure, achieving up to 70% compression. Likewise, Lindstrom’s ZFP compresses data in local blocks (Li14, ) with 1.3–2.6$\times$ (23–61%) compression. More recently, Di and Cappello (DiCa16, ) report up 3–436$\times$ (66–99.8%) compression, but again focusing on local blocks.

Our method, in contrast, aims at detecting global structure in the data, yielding much higher compression ratios in exchange for potential loss in local accuracy. It does not process the data in blocks but rather considers the data in its entirety. Before we delve into the mathematical details, we note briefly that lossy compression need not replace the original data; rather, the compressed version is a thumbnail or preview of the full dataset, which may reside on long-term storage or be regenerated.

In this paper, we consider the Tucker tensor decomposition (Tu66, ), also known as the higher-order singular value decomposition (HOSVD) (DeDeVa00, ). This is an effective tool for compression for many application domains (VaTe02, ; KaYaMe12, ; HaCaTe12, ; AfGiTa14, ; BaPa15, ; GaSoVe16, ). The idea is to consider the multiway structure of the data. For instance, the problem described above is a five-way object, and so we exploit this structure in the compression procedure. Specifically, each mode is compressed individually by determining a small number of vectors that span the fibers in that mode, fibers being the vectors that are the higher-order analogue of matrix rows and columns. We aim to compress an $N$-way tensor of size ${I}_{0}^{\phantom{{}^{\prime}}}\times{I}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{I}_{N-1}^{\phantom{{}^{\prime}}}$ to size ${R}_{0}^{\phantom{{}^{\prime}}}\times{R}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{R}_{N-1}^{\phantom{{}^{\prime}}}$ by computing $R_{n}$ vectors for each mode $n$ that approximately span the range of the mode-$n$ fibers. We call $R_{n}$ the rank of mode $n$. The ranks are typically selected to retain a specified relative accuracy $\epsilon$, i.e., if ${\bm{{\mathscr{{X}}}}}$ is the data tensor and ${\bm{\hat{\mathscr{{X}}}}}$ is the reconstruction from the compressed representation, then we choose the ranks such that

[TABLE]

Choosing $R_{n}=I_{n}$ yields perfect reconstruction, so viable choices for the ranks always exist. If we assume $I_{n}=I$ and $R_{n}=R$ for all $n$, for sake of exposition, then the storage is reduced from $I^{N}$ to $R^{N}+NIR$, and so the compression ratio is $\approx(I/R)^{N}$, which is the size of the original data divided by the size of the compressed representation.

In this paper, we develop a parallel implementation of the sequentially-truncated HOSVD (ST-HOSVD) (VaVaMe12, ). This papers builds on past work by Austin, Ballard, and Kolda (AuBaKo16, ), which showed initial results for parallel versions of ST-HOSVD as well as the higher-order orthogonal iteration (HOOI). Our contributions are as follows:

1. (1)

We describe the details of the TuckerMPI software which were not provided in (AuBaKo16, ), including global data distributions (even in the case when a tensor dimension is not divisible by the number of processors in that dimension), local data layouts, and sequential and parallel algorithms. 2. (2)

We present a new and improved kernel for computing the Gram matrix in ST-HOSVD, which is faster and less sensitive to the parallel data layout than the version in (AuBaKo16, ). 3. (3)

We present more extensive compression results for terabyte-scale data by compressing data sets that are an order of magnitude larger than those in (AuBaKo16, ). Specifically, we show that we can compress 4.5 and 6.7 TB datasets by 2–5 orders of magnitude in O(10-100) seconds, less time then reading the data from the parallel file system. 4. (4)

An advantage of the Tucker decomposition that was mentioned but not implemented for (AuBaKo16, ) is that we can reconstruct just a portion of the full tensor. Here, we present an efficient method for partial reconstruction, taking care not to create any object that is larger than the compressed or reconstructed subtensor so that it can be run on a workstation rather than a parallel system. We give experimental results to showcase its efficiency.

2. Notation and Mathematical Background

We use boldface Euler script letters to denote tensors (${\bm{{\mathscr{{X}}}}}$) and boldface uppercase letters to denote matrices (${\bm{\mathbf{{X}}}}$). We reserve the uppercase letters $I,J,K,L,M,N,P,R$ to denote sizes and the corresponding lowercase letters $i,j,k,l,m,n,p,r$ to denote the indices. We use zero-indexing throughout so that if $i$ is the index corresponding to size $I$, then we have $i=0,1,\dots,I-1$. For any size $I$, we use the notation $[\,I\,]$ to denote the set $\set{0,1,\dots,I-1}$.

2.1. Sizes

We define a few special quantities with respect to tuples of sizes, which are used in describing both tensor and processor grid sizes. For an object with dimensions ${I}_{0}^{\phantom{{}^{\prime}}}\times{I}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{I}_{N-1}^{\phantom{{}^{\prime}}}$, we define the product of all its sizes (the total size) as

[TABLE]

We further define some quantities that depend on the mode $n\in[\,N\,]$:

[TABLE]

For the edge cases, we say ${I}_{0}^{\varolessthan}=1$ and ${I}_{N-1}^{\varogreaterthan}=1$.

2.2. Tensor Operations

We discuss key tensor operations for computing the Tucker decomposition. Here we assume an $N$-way tensor ${\bm{{\mathscr{{X}}}}}$ of size ${I}_{0}^{\phantom{{}^{\prime}}}\times{I}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{I}_{N-1}^{\phantom{{}^{\prime}}}$.

2.2.1. Tensor Unfolding

The mode-$n$ unfolding of ${\bm{{\mathscr{{X}}}}}$ rearranges the elements of the $N$-way tensor into a matrix, denoted ${\bm{\mathbf{{X}}}}_{(n)}$, of size $I_{n}\times{I}_{n}^{\oslash}$. We map tensor element $(i_{0},i_{1},\dots,i_{N-1})$ to matrix element $(i_{n},{i}_{n}^{\prime})$ where

[TABLE]

See also Sec. 4.2 for examples of unfolded tensors and details of the organization in computer memory.

2.2.2. Tensor Norm

The norm of a tensor is the square root of the sum of squares of all the elements. This means it is equivalent to the Frobenious norm of any unfolding, i.e., $\|{\bm{{\mathscr{{X}}}}}\|=\|{\bm{\mathbf{{X}}}}_{(n)}\|_{F}$ for any $n\in[\,N\,]$.

2.2.3. Tensor Times Matrix

The mode-$n$ product of ${\bm{{\mathscr{{X}}}}}$ with a matrix ${\bm{\mathbf{{U}}}}$ of size $J\times I_{n}$ is denoted ${\bm{{\mathscr{{X}}}}}\times_{n}{\bm{\mathbf{{U}}}}$, and the result is of size $I_{0}\times\cdots\times I_{n-1}\times J\times I_{n+1}\times\cdots\times I_{N-1}$. This can be expressed in terms of unfolded tensors, i.e.,

[TABLE]

This is also known as the tensor-times-matrix (TTM) product. For different modes, the order is irrelevant so that

[TABLE]

In the same mode, order matters so that

[TABLE]

where ${\bm{\mathbf{{V}}}}$ is $K\times J$.

3. Review of the Sequentially-truncated higher-order SVD

In this section, we describe the Tucker decomposition and review the ST-HOSVD method that we parallelize to compute it. Before we do so, we introduce some notation and basic theory.

Let the $N$-way tensor ${\bm{{\mathscr{{X}}}}}$ of size ${I}_{0}^{\phantom{{}^{\prime}}}\times{I}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{I}_{N-1}^{\phantom{{}^{\prime}}}$ denote the data tensor to be compressed. The goal is to approximate ${\bm{{\mathscr{{X}}}}}$ as

[TABLE]

The tensor ${\bm{{\mathscr{{G}}}}}$ is called the core tensor, and its size is denoted by ${R}_{0}^{\phantom{{}^{\prime}}}\times{R}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{R}_{N-1}^{\phantom{{}^{\prime}}}$. Each factor matrix ${\mathbf{U}}_{n}$ is necessarily of size $I_{n}\times R_{n}$ for $n\in[\,N\,]$, and we assume throughout that the factor matrices have orthonormal columns. This means ${\mathbf{U}}_{n}^{\sf T}{\mathbf{U}}_{n}={\bm{\mathbf{{I}}}}$ (the $R_{n}\times R_{n}$ identity matrix). The storage of ${\bm{{\mathscr{{X}}}}}$ is $I^{\oast}$ as compared to the storage for ${\bm{\hat{\mathscr{{X}}}}}$ which is $R^{\oast}+\sum_{n\in[\,N\,]}R_{n}I_{n}$. If, for example, $I_{n}/R_{n}=2$ for all $n$, then the compression ratio is $I^{\oast}/R^{\oast}\approx 2^{N}$.

We review a few relevant facts about Tucker per (KoBa09, , Section 4.2). If the factor matrices are given, then it can be shown that the optimal ${\bm{{\mathscr{{G}}}}}$ is

[TABLE]

Substituting eq. 8 back into eq. 7 and using identity eq. 6, we have

[TABLE]

Thus, the factor matrices determine the projections (one per mode) of the original tensor down to the reduced space. It is easy to show that $\|{\bm{{\mathscr{{X}}}}}-{\bm{\hat{\mathscr{{X}}}}}\|^{2}=\|{\bm{{\mathscr{{X}}}}}\|^{2}-\|{\bm{{\mathscr{{G}}}}}\|^{2}$, which means that ${\bm{{\mathscr{{G}}}}}$ retains most of the mass of ${\bm{{\mathscr{{X}}}}}$ but is just represented with respect to a different basis.

It is convenient to make a few special definitions. Assume that the factor matrices are specified. Then define the partial core that is the result of applying $n+1$ factor matrices to ${\bm{{\mathscr{{X}}}}}$:

[TABLE]

This means that ${\bm{\tilde{\mathscr{{G}}}}}_{n}$ is only reduced in the first $(n+1)$ modes and so is of size $R_{0}\times\cdots\times R_{n}\times I_{n+1}\cdots\times I_{N-1}$. By definition, ${\bm{{\mathscr{{G}}}}}={\bm{\tilde{\mathscr{{G}}}}}_{N-1}$. Likewise, we can define the incremental approximation using the partial core to be

[TABLE]

By definition, ${\bm{\hat{\mathscr{{X}}}}}={\bm{\tilde{\mathscr{{X}}}}}_{N-1}$.

3.1. Key Theorem

The following theorem categorizes the error in terms of the incremental approximations and orthogonal projections based on the factor matrices. We refer the reader to (VaVaMe12, ) for the proof. This theorem is key to understanding how to pick the factor matrices in the ST-HOSVD.

Theorem 3.1 (Vannieuwenhoven, Vandebril, and Meerbergen (VaVaMe12, , Theorem 5.1)).

Let ${\bm{{\mathscr{{X}}}}}$ be a tensor of size ${I}_{0}^{\phantom{{}^{\prime}}}\times{I}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{I}_{N-1}^{\phantom{{}^{\prime}}}$ which is approximated by ${\bm{\hat{\mathscr{{X}}}}}$ as defined in eq. 9 in which the factor matrices ${\mathbf{U}}_{n}$ of size $R_{n}\times I_{n}$ have orthonormal columns. The approximation error is then given by

[TABLE]

where ${\bm{\tilde{\mathscr{{X}}}}}_{n}$ is as defined in eq. 11. Additionally, the error is bounded by

[TABLE]

3.2. HOSVD Method

The bound in eq. 13 from Theorem 3.1 is key to understanding the HOSVD method, originally known as Tucker1 (Tu66, ; DeDeVa00, ). Suppose the goal is to find a Tucker decomposition of the form in eq. 7 with relative error no greater then $\epsilon$, i.e.,

[TABLE]

Then the idea is as follows. Let the eigenvalue decomposition of the Gram matrix ${\bm{\mathbf{{S}}}}$ of the mode-$n$ unfolding be given by

[TABLE]

Here, ${\bm{\mathbf{{\Lambda}}}}=\operatorname{diag}(\set{\lambda_{1},\dots,\lambda_{I_{n}}})$ and $\lambda_{1}\geq\lambda_{2}\cdots\geq\lambda_{I_{n}}\geq 0$ are the eigenvalues in descending order. The matrix ${\bm{\mathbf{{V}}}}$ contains the corresponding eigenvectors. We choose ${\mathbf{U}}_{n}$ and $R_{n}$ so that

[TABLE]

This choice of ${\mathbf{U}}_{n}$ ensures that

[TABLE]

Repeating this procedure for each $n\in[\,N\,]$ ensures that the desired error bound eq. 14 holds. This discussion can be framed equivalently in terms of the leading left singular vectors of ${\bm{\mathbf{{X}}}}_{(n)}$, which is how the HOSVD is usually described. However, we explicitly form the Gram matrices in our parallelized method, so this presentation is more convenient in our exposition.

3.3. ST-HOSVD

The HOSVD works with the full tensor at each step. The idea behind the sequentially-truncated version introduced in (VaVaMe12, ) is that we work with the partial cores. In other words, at step $n$, we compute the next factor matrix based on ${\bm{\tilde{\mathscr{{G}}}}}_{n-1}$. This is possible because we can swap ${\bm{\tilde{\mathscr{{G}}}}}_{n-1}$ for ${\bm{\tilde{\mathscr{{X}}}}}_{n-1}$ in the summands in the error expression eq. 12 due to the following equivalence:

[TABLE]

Hence, substituting eq. 17 into eq. 12 then yields the following key corollary.

Corollary 3.2.

Let the conditions of Theorem 3.1 hold. Then

[TABLE]

where ${\bm{\tilde{\mathscr{{G}}}}}_{n}$ is as defined in eq. 10.

Corollary 3.2 means that we can pick the $(n+1)$st factor matrix, ${\mathbf{U}}_{n}$, using the partial core ${\bm{\tilde{\mathscr{{G}}}}}_{n-1}$. Choosing ${\mathbf{U}}_{0}$ is unchanged. For $n>1$, however, we replace the eigenvalue problem in eq. 15 with

[TABLE]

We can still pick ${\mathbf{U}}_{n}$ according to eq. 16, but we are just working with the eigen-decomposition of a different matrix. The ST-HOSVD algorithm is given in Alg. 1.

3.4. Quasi-Optimality

Neither the HOSVD nor the ST-HOSVD yields an optimal rank-$(R_{0},R_{1},\dots,R_{N-1})$ decomposition, and neither is necessarily more accurate than the other. The major advantage of the ST-HOSVD is in terms of computational cost: the cost of the Gram computation is decreased by a factor of $I_{n}/R_{n}$ at each step. Both methods yield quasi-optimal results with error within a factor of the square root of the number of modes of the best approximation, as specified by the following theorem from (Hackbusch14, ).

Theorem 3.3 (Hackbusch (Hackbusch14, , Theorems 6.9-10)).

Let ${\bm{{\mathscr{{X}}}}}$ be a tensor of size ${I}_{0}^{\phantom{{}^{\prime}}}\times{I}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{I}_{N-1}^{\phantom{{}^{\prime}}}$, ${\bm{\hat{\mathscr{{X}}}}}_{\rm{HOSVD}}$ be the rank ${R}_{0}^{\phantom{{}^{\prime}}}\times{R}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{R}_{N-1}^{\phantom{{}^{\prime}}}$ approximation computed by HOSVD, and ${\bm{\hat{\mathscr{{X}}}}}_{\rm{ST-HOSVD}}$ be the rank ${R}_{0}^{\phantom{{}^{\prime}}}\times{R}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{R}_{N-1}^{\phantom{{}^{\prime}}}$ approximation computed by ST-HOSVD. Then

[TABLE]

and

[TABLE]

where ${\bm{\hat{\mathscr{{X}}}}}_{\rm{opt}}$ is the optimal rank-${R}_{0}^{\phantom{{}^{\prime}}}\times{R}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{R}_{N-1}^{\phantom{{}^{\prime}}}$ approximation of ${\bm{{\mathscr{{X}}}}}$.

3.5. Pre- and Post-Processing

Algorithms for computing Tucker models like ST-HOSVD (Alg. 1) minimize the relative approximation error in a norm-wise sense. As a result, the relative component-wise errors are generally smaller for large entries in the data tensor than for small entries. If the data is not preprocessed, the differences in magnitude are artifacts of the type of variable or unit of measurement rather than reflecting the importance of the data value.

To address these issues, we pre-process the original data using hyperslice-wise computations, where a single hyperslice corresponds to the values of a particular variable like pressure or temperature, all with the same physical units. These computations include gathering statistics, such as mean or maximum value, for each hyperslice of a particular mode and applying a single linear function to every element in a hyperslice to uniformly shift and/or scale the values. The vector of shifting and scaling values are stored so that the inverse operations can be applied in post-processing after the (partial) reconstruction process.

The hyperslice-wise statistics that can be collected include the mean, standard deviation, maximum, minimum, vector 1-norm, and vector 2-norm. Common preprocessing techniques are (a) shifting by the mean value (to center each hyperslice at 0) and scaling by the standard deviation (to impose a variance of

1. and (b) scaling by the maximum absolute value (to impose a range of ${-}1$ to 1). We refer to (a) as “standardization” and (b) as “max rescaling.”

4. Data Layouts

In this section, we describe how tensors are stored in local memory and how they are distributed for parallel computation. Factor matrices are stored redundantly on every processor. We assume throughout that we have a generic $N$-way tensor ${\bm{{\mathscr{{Y}}}}}$ of size ${J}_{0}^{\phantom{{}^{\prime}}}\times{J}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{J}_{N-1}^{\phantom{{}^{\prime}}}$. WLOG, this ${\bm{{\mathscr{{Y}}}}}$ tensor is the one that is updated in Alg. 1, so each $J_{n}$ is either $R_{n}$ or $I_{n}$. We assume that the processors are logically arranged into an $N$-way Cartesian processor grid, with dimensions ${P}_{0}^{\phantom{{}^{\prime}}}\times{P}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{P}_{N-1}^{\phantom{{}^{\prime}}}$, as in (AuBaKo16, ; SK16, ; BHK18-TR, ). In the analysis, we assume $P_{n}\in\set{1,\dots,J_{n}}$ for each $n\in[\,N\,]$. We use the size shorthands discussed in Sec. 2.1, e.g., $P^{\oast}$ denotes the total number of processors.

4.1. Tensor Data Layout

Consider the natural descending format that generalizes the column-major matrix layout. That is, tensor element $(j_{0},j_{1},\dots,j_{N-1})$ can be mapped to

[TABLE]

which we refer to as the linear index. We can also define the inverse operation,

[TABLE]

Fig. 1 shows a $3\times 4\times 3\times 2$ tensor with each entry labeled by its linear index.

4.2. Unfolded Tensor Data Layout

If the tensor is stored in the natural descending format, then the data layout of each unfolding can be thought of as a set of contiguous submatrices, where each submatrix is stored in row-major ordering. In other words, there is no reason to do any data movement in memory to work with the unfolded matrix because BLAS calls can operate on submatrices stored in row- or column-major order, as observed in (AuBaKo16, ; LBPSV15, ; HBJT18, ; PTC13a, ). The number and dimensions of these submatrices depend on the mode of the unfolding. For the $n$th-mode unfolding, the number of submatrices is ${J}_{n}^{\varogreaterthan}$, and each submatrix is of size ${J}_{n}^{\phantom{{}^{\prime}}}\times{J}_{n}^{\varolessthan}$. Let’s look at these for the $3\times 4\times 3\times 2$ tensor in Fig. 1. Recall that the values indicate the relative position in the natural descending format. In the case $n=0$, each submatrix has one column, which corresponds to a global column-major matrix, so it has 24 submatrices of size $3\times 1$:

[TABLE]

Note that this yields a special mode-0 case in the implementation of operations because it is more efficient to treat the data as one column-major block than many blocks of vectors. The mode-1 unfolding has 6 row-major submatrices of size $4\times 3$:

[TABLE]

The mode-2 unfolding has 2 row-major submatrices of size $3\times 12$:

[TABLE]

In the case $n=N-1$, there is only one submatrix, which corresponds to a global row-major matrix, so we have one submatrix of size $2\times 36$:

[TABLE]

These layouts become important when we discuss the local operations in Sec. 5.

4.3. Global Tensor Data Distribution

Now we consider how the tensor is distributed on the global processor grid. We let an overbar indicate local quantities with respect to a specific processor. If $J_{n}$ is the global size of mode $n$, then $\bar{J}_{n}$ is the corresponding local size for the subtensor owned by a processor $(\bar{p}_{0},\bar{p}_{1},\dots,\bar{p}_{N-1})$.

Processor $(\bar{p}_{0},\bar{p}_{1},\dots,\bar{p}_{N-1})$ owns a block of size $\bar{J}_{0}^{\phantom{{}^{\prime}}}\times\bar{J}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times\bar{J}_{N-1}^{\phantom{{}^{\prime}}}$ where $\bar{J}_{n}$ is given by

[TABLE]

Tab. 1 shows the size of each local subtensor for a tensor of size $3\times 4\times 3\times 2$ distributed on a processor grid of size $2\times 2\times 2\times 1$. Tab. 1 also indicates the processor linear index, which is defined analogously to eq. 20:

[TABLE]

The tensor element with global index $(j_{0},j_{1},\dots,j_{N-1})$ is mapped to processor $(\bar{p}_{0},\bar{p}_{1},\dots,\bar{p}_{N-1})$ where $\bar{p}_{n}$ is given by

[TABLE]

Fig. 2a shows the mapping of elements to processors for the tensor in Fig. 1 using a processor grid of size $2\times 2\times 2\times 1$.

We define the set of mode-$n$ global indices mapped to processor $\bar{p}_{n}$ to be

[TABLE]

The per-mode assignments can be combined to get all the indices assigned to a specific processor. One interesting note is that elements that are contiguous in the global natural descending linear index are not necessarily mapped to the same processor. For instance, in the example in Fig. 2, consider that processor $\bar{p}^{\prime}=2$ owns the following global linear indices:

[TABLE]

which is not contiguous in the global linear ordering.

The local index is denoted $(\bar{j}_{0},\bar{j}_{1},\dots,\bar{j}_{N-1})$ where $j_{n}$ is given by

[TABLE]

Locally, the subtensors are stored contiguously according to the natural descending index and in the same order as their global linear index. Fig. 2b shows the local linear indices for the example in Fig. 2a. This equation can be reversed to find the global index $(j_{0},j_{1},\dots,j_{N-1})$ where $j_{n}$ is given by

[TABLE]

4.4. Global Unfolding Data Distribution

In the parallel case, the unfolded tensor is distributed among the processors in a block fashion. For instance, the mode-1 unfolding is here color-coded to the processor (using the same colors as in Fig. 2):

[TABLE]

Note that a single MPI process may own multiple contiguous pieces of the unfolded tensor. In particular, the data distribution of the unfoldings are not a standard distribution such as blocked or block-cyclic. However, the distribution is still a blocked matrix distribution, equivalent to unfolding the local tensor. Here we show the local linear index of each entry:

[TABLE]

Notice that they are contiguous for each processor. This means that we can work with the locally unfolded tensor for the local operations in distributed computations.

4.5. Processor Fibers

Just as for the data tensor, there is a corresponding processor grid “unfolding” that corresponds to the tensor unfoldings. Instead of an $N$-way processor grid, we can think of the processors as rearranged into a 2-way grid. So, if we revisit the example in the previous subsection, the mode-1 processor unfolding for the $2\times 2\times 2\times 1$ processor grid is:

[TABLE]

We see the processors in the same column match in every index except the mode-1 index, since this is the mode-1 unfolding. We see later that certain operations require collective communications within each processor column in this unfolding. We refer to these column groups of processors as mode-$n$ processor fibers. We refer to row groups of processors as mode-$n$ processor slices.

5. Local Kernels

In Alg. 1, there are three key kernels: Gram of the mode-$n$ unfolded tensors in 4, the $I_{n}\times I_{n}$ eigenproblem in 7, and tensor-times-matrix to shrink in mode $n$ in 8. Here we explain the serial implementation, which is also used for the local computations in the parallelized version. We explain all the functions with respect to a generic tensor ${\bm{{\mathscr{{Y}}}}}$ of size ${J}_{0}^{\phantom{{}^{\prime}}}\times{J}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{J}_{N-1}^{\phantom{{}^{\prime}}}$.

5.1. Local Gram

We want to compute ${\bm{\mathbf{{S}}}}={\bm{\mathbf{{Y}}}}_{(n)}^{\phantom{T}}{\bm{\mathbf{{Y}}}}_{(n)}^{\sf T}$ so ${\bm{\mathbf{{S}}}}$ will be of size $J_{n}\times J_{n}$. The arithmetic cost is $\mathcal{O}(J^{\oast}J_{n})$. Although the computation is mathematically one matrix operation, the data layout of the unfolding of the input tensor prevents a single call to the syrk BLAS subroutine, which requires strided row- or column-major ordering. Thus, the algorithm works on the native row-major submatrices as described in Sec. 4.2, except in the case of $n=0$ where the algorithm works on the native column-major matrices. For $n>0$, ${\bm{\mathbf{{Y}}}}_{(n)}$ has ${J}_{n}^{\varogreaterthan}$ submatrices of size $J_{n}\times{J}_{n}^{\varolessthan}$. We denote the $j$th block-column submatrix as ${\bm{\mathbf{{Y}}}}_{(n)}\left[j\right]\equiv{\bm{\mathbf{{Y}}}}_{(n)}(\,:\,,j\!\cdot\!{J}_{n}^{\varolessthan}:j\!\cdot\!{J}_{n}^{\varolessthan}+{J}_{n}^{\varolessthan}-1)$, which is stored in row-major form. The algorithm is shown in Alg. 2.

5.2. Local Eigenvalue Decomposition

We need to compute the eigenvalue decomposition of a $J_{n}\times J_{n}$ matrix ${\bm{\mathbf{{S}}}}$. The computational cost is $\mathcal{O}(J_{n}^{3})$. We use the the LAPACK direct eigenvalue computation routine syev to compute all the eigenvalues and all the eigenvectors.

A couple of notes are in order. The cost of the full eigenvector decomposition is approximately $\frac{10}{3}J_{n}^{3}$. One alternative approach would be to compute the full set of eigenvalues at a cost of $\frac{4}{3}J_{n}^{3}$, and then compute only the leading eigenvectors with $O(J_{n}^{2}R_{n})$ flops. Iterative methods, such as subspace iteration, would also work well in this case. However, because this phase of computation has never been a bottleneck for our applications, we have not implemented these cheaper approaches. If $J_{n}$ is relatively large compared to the product of the other dimensions, ${J}_{n}^{\oslash}$, then the eigenvalue decomposition may become a bottleneck, so we leave this as a topic for future work.

5.3. Local TTM

We want to compute the product ${\bm{{\mathscr{{Z}}}}}={\bm{{\mathscr{{Y}}}}}\times_{n}{\bm{\mathbf{{V}}}}$ where ${\bm{\mathbf{{V}}}}$ is a matrix of size $K\times J_{n}$. Recall that the TTM operation is defined in Sec. 2.2.3. The arithmetic cost of TTM is $\mathcal{O}(J^{\oast}K)$.

As in the case of Gram, although TTM is mathematically a matrix multiplication (${\bm{\mathbf{{Z}}}}_{(n)}={\bm{\mathbf{{V}}}}{\bm{\mathbf{{Y}}}}_{(n)}$), the data layouts of ${\bm{\mathbf{{Y}}}}_{(n)}$ and ${\bm{\mathbf{{Z}}}}_{(n)}$ prevent a single call to the gemm subroutine in BLAS because BLAS requires strided row- or column-major access to the matrices. We use the same notation as for Gram, so that ${\bm{\mathbf{{Y}}}}_{(n)}\left[j\right]$ is the $j$th block column of ${\bm{\mathbf{{Y}}}}_{(n)}$ stored in row-major form. Likewise, ${\bm{\mathbf{{Z}}}}_{(n)}$ is organized into row-major block column submatrices of size $K\times{J}_{n}^{\varolessthan}$, and the $j$th submatrix is denoted as ${\bm{\mathbf{{Z}}}}_{(n)}\left[j\right]$. In the case of $n=0$, both ${\bm{\mathbf{{Y}}}}_{(n)}$ and ${\bm{\mathbf{{Z}}}}_{(n)}$ are natively in column-major mode. The algorithm is shown in Alg. 3.

6. Distributed Kernels

To parallelize STHOSVD (Alg. 1), we use parallel algorithms for the two key kernels: Gram and TTM. All other computations are performed redundantly on each processor, including the eigendecomposition.

Throughout this section, we consider a generic tensor ${\bm{{\mathscr{{Y}}}}}$ of size ${J}_{0}^{\phantom{{}^{\prime}}}\times{J}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{J}_{N-1}^{\phantom{{}^{\prime}}}$, distributed on a processor grid of size ${P}_{0}^{\phantom{{}^{\prime}}}\times{P}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{P}_{N-1}^{\phantom{{}^{\prime}}}$, as described in Sec. 4.3. We use an overbar to denote the local portions/versions/sizes of distributed variables. For instance, ${\bm{\bar{\mathscr{{Y}}}}}$ denotes the local portion of ${\bm{{\mathscr{{Y}}}}}$, and $\bar{J}_{0}^{\phantom{{}^{\prime}}}\times\bar{J}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times\bar{J}_{N-1}^{\phantom{{}^{\prime}}}$ is its size. Note that local quantities, including their sizes, may vary from processor to processor.

6.1. Assumptions on Collective Communication

To analyze our algorithms, we use the MPI model of distributed-memory parallel computation. We assume a fully connected network of $P$ processors and therefore do not model network contention. For simplicity of discussion in our analysis, we assume that $P$ is a power of two and that optimal collective communication algorithms are used. The cost to send/receive a message of size $W$ words between any two processors is $\alpha+W\beta$, where $\alpha$ is the latency cost and $\beta$ is the per-word transfer cost. The cost of the collective communications used in this work are given in Tab. 2, with $\gamma$ corresponding to the time per floating point operation (flop). For simplicity of presentation, we will ignore the flop cost of the reductions in later analysis, as they are typically dominated by the bandwidth costs. For more discussion of the model and descriptions of efficient collectives, see (CH+07, ; TRG05, ).

6.2. Parallel Gram

The goal here is to compute ${\bm{\mathbf{{S}}}}={\bm{\mathbf{{Y}}}}_{(n)}^{\phantom{T}}{\bm{\mathbf{{Y}}}}_{(n)}^{\sf T}$ where ${\bm{{\mathscr{{Y}}}}}$ is distributed as described in Sec. 4.3. Each processor owns ${\bm{\bar{\mathscr{{Y}}}}}$ (local portion of ${\bm{{\mathscr{{Y}}}}}$) of size $\bar{J}_{0}^{\phantom{{}^{\prime}}}\times\bar{J}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times\bar{J}_{N-1}^{\phantom{{}^{\prime}}}$ where $\bar{J}_{n}=\operatorname{lsz}\{\,\bar{p}_{n},J_{n},P_{n}\}$ and the mode-$n$ indices correspond to the global indices in $\bar{\mathcal{J}}_{n}=\operatorname{prcmap}\{\,\bar{p}_{n},{J_{n}},P_{n}\}$. In the end, each processor will own the entirety of ${\bm{\mathbf{{S}}}}$, which is only of size $J_{n}\times J_{n}$.

Austin et al. (AuBaKo16, ) previously proposed a parallel Gram as illustrated in Fig. 3. In this algorithm, each processor fiber, which owns a column block of ${\bm{\mathbf{{Y}}}}_{(n)}$, computes a contribution ${\bm{\mathbf{{V}}}}$ to the result ${\bm{\mathbf{{S}}}}$; ${\bm{\mathbf{{V}}}}$ is distributed across the processor fiber. In order to compute ${\bm{\mathbf{{V}}}}$ within the fiber, the processors rotate their tensor data around in a round-robin fashion, and at each step, each processor computes a $\bar{J}_{n}^{\phantom{{}^{\prime}}}\times\bar{J}_{n}^{\phantom{{}^{\prime}}}$ block of ${\bm{\mathbf{{V}}}}$. In order to compute ${\bm{\mathbf{{S}}}}$, the processors perform an All-Reduce across processor slices, so that ${\bm{\mathbf{{S}}}}$ is redundantly stored on every processor fiber but distributed across the processors within the fiber. We refer to the older version as the round-robin variant.

We propose a new version of parallel Gram that is nearly always faster, which we refer to as the redistribution variant. Alg. 4 presents the new parallel algorithm for computing the Gram matrix corresponding to a particular mode. The algorithm assumes that the input tensor ${\bm{{\mathscr{{Y}}}}}$ is block distributed; at the end of the algorithm, the output matrix ${\bm{\mathbf{{S}}}}={\bm{\mathbf{{Y}}}}_{(n)}^{\phantom{T}}{\bm{\mathbf{{Y}}}}_{(n)}^{\sf T}$ is redundantly stored on every processor. The parallel Gram algorithm presented here differs from the one described in (AuBaKo16, ), depicted in Fig. 3. When the number of processors in the specified mode $P_{n}$ is 1, the algorithms are identical. However, when $P_{n}>1$, the previous algorithm uses $P_{n}-1$ communication steps within the processor fiber and then communicates across the processor slice. As we describe in more detail below, Alg. 4 works by first redistributing the tensor data with an All-to-All collective within the processor fiber and then performing a reduction across all processors. We compare the communication costs of the two algorithms at the end of the section.

As in the case of Alg. 5, 2 and 3 of Alg. 4 define the processor’s index and the set of processors within the processor’s $n$th-mode processor fiber. The goal of 5 is to obtain a 1D parallel distribution of the tensor. With this distribution, all processors can perform Gram computations with their local data, and the only remaining communication is to sum up the results over all processors. Fig. 4 illustrates the redistribution. Because each processor fiber stores a set of columns of the matricized tensor (distributed row-wise), the redistribution occurs within each fiber independently and converts the row-wise distribution to a column-wise distribution.

Again, a key issue in the implementation is the need to pack and unpack buffers for the All-to-All collective, depending on $n$. The input buffer must be arranged so that the data to be received by each processor is stored contiguously, and the result buffer is ordered so that the data each processor sent is contiguous. In the case $n=0$, the local matricized tensor is column-major, so no packing is necessary. After the All-to-All, the result buffer consists of ${P}_{n}^{\phantom{{}^{\prime}}}$ contiguous column-major blocks, where each column has length $\bar{J}_{n}^{\phantom{{}^{\prime}}}$. The unpacking consists of collecting the ${P}_{n}^{\phantom{{}^{\prime}}}$ chunks of each mode-$n$ fiber (of length ${J}_{n}^{\phantom{{}^{\prime}}}$) into contiguous columns. In the case $0<n<N-1$, the local matricized tensor is stored as an array of row-major submatrices. The packing consists of converting every row-major submatrix to column-major ordering, which makes every local mode-$n$ fiber contiguous. The unpacking is the same as for $n=0$: for each mode-$n$ fiber, the ${P}_{n}^{\phantom{{}^{\prime}}}$ chunks are made contiguous. We note that the local ordering of the columns is not consistent with a matricized tensor format, but the Gram computation is invariant under column permutations, so this does not affect the result. In the case $n=N-1$, the local matricized tensor is row-major. Instead of converting the entire matrix to column-major, each contiguous row is broken up into ${P}_{n}^{\phantom{{}^{\prime}}}$ chunks, and the row ordering is maintained. After the All-to-All, no unpacking is necessary because the result buffer is a (vertical) concatenation of row-major matrices and is thus also row-major.

After obtaining a 1D distribution of the matricized tensor, each processor performs its local computation in 6, computing the Gram matrix associated with its subset of the columns of the matricized tensor. In order to compute the Gram matrix of the entire matricized tensor, all processors participate in an All-Reduce (7), which results in (symmetric) ${\bm{\mathbf{{S}}}}$ being redundantly stored on all processors.

The cost of Alg. 4 is

[TABLE]

For comparison, the cost of the previous parallel Gram algorithm (AuBaKo16, ) is

[TABLE]

where the major difference is in the bandwidth cost. In particular, Alg. 4 is more efficient when $P_{n}\bar{J}^{\oast}\geq J_{n}^{2}$. In this case, both bandwidth cost terms of the new algorithm are smaller than the first term of the old algorithm. (If the converse is true, then both bandwidth cost terms of the old algorithm are smaller than the second term of the new algorithm.) We expect the inequality to hold in nearly all cases (except for extreme strong-scaling cases) because $\bar{J}^{\oast}$ itself is the size of the local input tensor, while $J_{n}^{2}$ is the size of the $n$th-mode Gram matrix.

The temporary memory requirement of Alg. 4 (besides the input and output data) is $2\bar{J}^{\oast}+J_{n}^{2}$ words. Two temporary arrays of the size of the local tensor are required for the send and receive buffers in the All-to-All (local data has to be reordered to match the requirements of the input buffer), and $J_{n}^{2}$ space is required for ${\bm{\mathbf{{W}}}}$ (All-Reduce requires separate input and output buffers). For comparison, the temporary memory requirement of the previous Gram algorithm is $\bar{J}^{\oast}+\bar{J}_{n}^{2}$.

6.3. Parallel TTM

We consider the problem of computing the TTM ${\bm{{\mathscr{{Z}}}}}={\bm{{\mathscr{{Y}}}}}\times_{n}{\bm{\mathbf{{V}}}}$ where the input ${\bm{{\mathscr{{Y}}}}}$ of size ${J}_{0}^{\phantom{{}^{\prime}}}\times{J}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{J}_{N-1}^{\phantom{{}^{\prime}}}$ is block distributed, the matrix ${\bm{\mathbf{{V}}}}$ of size $K_{n}\times J_{n}$ with $K_{n}<J_{n}$ is stored redundantly on every processor, and the output tensor ${\bm{{\mathscr{{Z}}}}}$ of size $J_{1}\times\cdots\times J_{n-1}\times K_{n}\times J_{n+1}\times\cdots\times J_{N}$ will be block distributed. Specifically, the data is distributed so that processor $(\bar{p}_{0},\bar{p}_{1},\dots,\bar{p}_{N-1})$ owns the following:

• •

${\bm{\bar{\mathscr{{Y}}}}}$ (local portion of ${\bm{{\mathscr{{Y}}}}}$) is of size $\bar{J}_{0}^{\phantom{{}^{\prime}}}\times\bar{J}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times\bar{J}_{N-1}^{\phantom{{}^{\prime}}}$ where $\bar{J}_{n}=\operatorname{lsz}\{\,\bar{p}_{n},J_{n},P_{n}\}$ and the mode-$n$ indices correspond to the global indices in $\bar{\mathcal{J}}_{n}=\operatorname{prcmap}\{\,\bar{p}_{n},{J_{n}},P_{n}\}$.

• •

${\bm{\bar{\mathscr{{Z}}}}}$ (local portion of ${\bm{{\mathscr{{Z}}}}}$) is of size $\bar{J}_{1}\times\cdots\times\bar{J}_{n-1}\times\bar{K}_{n}\times\bar{J}_{n+1}\times\cdots\times\bar{J}_{N}$ where $\bar{K}_{n}=\operatorname{lsz}\{\,\bar{p}_{n},K_{n},P_{n}\}$; it is distributed the same as ${\bm{\bar{\mathscr{{Y}}}}}$ except that the mode-$n$ indices correspond to the global indices in $\bar{\mathcal{K}}_{n}=\operatorname{prcmap}\{\,\bar{p}_{n},{K_{n}},P_{n}\}$.

• •

${\bm{\bar{\mathbf{{V}}}}}$ is of size $K_{n}\times\bar{J}_{n}$ and is the submatrix of ${\bm{\mathbf{{V}}}}$ corresponding to the columns in the set $\bar{\mathcal{J}}_{n}$. (Although every processor owns all of ${\bm{\mathbf{{V}}}}$, the distributed TTM only needs ${\bm{\bar{\mathbf{{V}}}}}$.)

Alg. 5 presents the parallel algorithm for distributed TTM. This is the same algorithm as previously presented by Austin et al. (AuBaKo16, ), but here we provide additional implementation details and analysis.

The computation reduces to a large matrix-matrix product, i.e., ${\bm{\mathbf{{Z}}}}_{(n)}={\bm{\mathbf{{V}}}}{\bm{\mathbf{{Y}}}}_{(n)}$. If ${\bm{\mathbf{{Y}}}}_{(n)}$ is partitioned into column blocks, the computation can be computed separately in each block. Therefore, since each mode-$n$ processor fiber owns a separate column block of ${\bm{\mathbf{{Y}}}}_{(n)}$, they are independent. For this reason, 3 defines the local processor fiber, and all communication for that processor occurs within the $P_{n}$ nodes of that fiber. There are ${P}_{n}^{\oslash}$ independent column fibers.

The algorithm chooses between two methods based on the size of $K_{n}$. The decision is based on the size of the intermediate quantities that are computed. The choice of method ensures that the temporary memory never exceeds the memory of ${\bm{\bar{\mathscr{{Y}}}}}$.

In the Reduce-Scatter variant, each processor computes the local matrix-matrix product, ${\bm{\bar{\mathbf{{W}}}}}={\bm{\bar{\mathbf{{V}}}}}{\bm{\bar{\mathbf{{Y}}}}}_{(n)}$, and then sums and distributes the result using a Reduce-Scatter. The temporary object ${\bm{\bar{\mathbf{{W}}}}}$ is of size $K_{n}\times\bar{J}_{n}^{\oslash}$, which may be too large to store on the processor. We assume that the processor has enough extra memory to store something the same size as ${\bm{\bar{\mathbf{{Y}}}}}_{(n)}$ (i.e., ${\bm{\bar{\mathscr{{Y}}}}}$), which is of size $\bar{J}_{n}\times{\bar{J}}_{n}^{\oslash}$, and so we use this variant if $K_{n}\leq\lfloor J_{n}/P_{n}\rfloor\approx\bar{J}_{n}$. This process is demonstrated in Fig. 5 where five processors are shaded (in the first processor fiber) to understand that data ownership of each processor.

A key issue in the implementation is that the entries’ ordering in ${\bm{\bar{\mathscr{{W}}}}}$ is not correct for a Reduce-Scatter and so must be reorganized when the data is packed for the Reduce-Scatter to obtain the proper layout of ${\bm{\bar{\mathscr{{Z}}}}}$ at the end of the call. The input buffer for Reduce-Scatter must be arranged so that contributions to each processor’s result ${\bm{\bar{\mathscr{{Z}}}}}$ are contiguous. The ordering of ${\bm{\bar{\mathscr{{W}}}}}$ consists of $\bar{J}_{n}^{\varogreaterthan}$ row-major submatrices of dimension $K\times\bar{J}_{n}^{\varolessthan}$, and it must be reordered into ${P}_{n}^{\phantom{{}^{\prime}}}$ contiguous subblocks, each consisting of $\bar{J}_{n}^{\varogreaterthan}$ row-major submatrices of dimension $K/{P}_{n}^{\phantom{{}^{\prime}}}\times\bar{J}_{n}^{\varolessthan}$ (assuming ${P}_{n}^{\phantom{{}^{\prime}}}$ divides $K$ evenly). If $n=N-1$ (the last dimension), no reordering is necessary, and no unpacking is necessary after the Reduce-Scatter for any dimension.

In the multiple-reduction variant, the algorithm uses a blocked approach that involves $P_{n}$ local TTMs and $P_{n}$ collective communications. Here, each iteration computes the contribution to one processor’s output (within the processor fiber — all processor fibers are working concurrently) and uses a Reduce collective to compute the sum across all processors in the fiber and store the result on the $\ell$th processor. The matrix ${\bm{\bar{\mathbf{{V}}}}}$ is divided into block rows so that ${\bm{\bar{\mathbf{{V}}}}}[\ell]$ owns the rows in the set $\operatorname{prcmap}\{\,\bar{p}_{n},{\ell,K_{n},P_{n}},P_{n}\}$. In this case, the dimensions of the temporary tensor ${\bm{\bar{\mathscr{{W}}}}}$ in the $\ell$th iteration is bounded above by $\bar{J}_{0}\times\cdots\times\lceil K_{n}/P_{n}\rceil\times\cdots\times\bar{J}_{N-1}$, which is essentially the same size as ${\bm{\bar{\mathscr{{Z}}}}}$ (recall that $\hat{K}_{\ell}$ is within 1 of $\bar{K}_{n}$). This process is demonstrated in Fig. 6, where we highlight the contribution to the 4th processor in the 1st column fiber.

The two variants perform the same number of flops and communicate almost the same amount of data (to within a factor of 2). The number of flops is $\mathcal{O}(J^{\oast}K_{n}/P^{\oast})$, which is the cost of the local TTM(s).

The data communication costs differ by only a factor of 2. In the Reduce-Scatter variant, the amount of data communicated is the cost of one Reduce-Scatter collective over the mode-$n$ fiber: $\prod_{m\neq n}\lceil J_{m}/P_{m}\rceil\lceil K_{n}/P_{n}\rceil(P_{n}-1)\approx{J}_{n}^{\oslash}K_{n}(P_{n}-1)/P^{\oast}$. In the multiple-reduction variant, there are $P_{n}$ Reduce collectives of a data size that is smaller by a factor of $P_{n}$, which yields the same asymptotic cost, but incurs an extra factor of 2 due to the cost of the Reduce collective.

The number of messages is fewer for the Reduce-Scatter variant, which has only one collective, at a cost of $\log P_{n}$ messages. The multiple-reduction variant involves $P_{n}$ collectives, at a total cost of $2P_{n}\log P_{n}$ messages.

The temporary memory (${\bm{\bar{\mathscr{{W}}}}}$) of the multiple-reduction variant is much lower: $\lceil K_{n}/P_{n}\rceil\bar{J}_{n}^{\oslash}$ words for the multiple-reduction variant versus $K_{n}\bar{J}_{n}^{\oslash}$ words for the Reduce-Scatter variant. As an aside, we note the block size of the blocked algorithm can be chosen arbitrarily, navigating a tradeoff between latency cost and memory footprint; we used a version that corresponded to the result size for simplicity, but it assumes the available remaining memory is at least the size of the required memory for the problem.

7. Parallel ST-HOSVD Cost Analysis

In this section we analyze the computation, communication, and temporary memory requirements of the ST-HOSVD algorithm. We note that the computation and bandwidth costs are sensitive to mode order, while the latency cost and memory requirements are not. This analysis assumes that the mode order used by the algorithm is increasing by mode index; the costs for other mode orders can be derived by relabeling modes. If the core ranks are specified a priori, then an optimal (in terms of flops or communication) mode ordering can be determined similar to the case of reconstruction, as described in Sec. 8. We also note that the communication costs are sensitive to the processor grid, but the computation cost and memory requirements are not. In this analysis, we use the new Gram algorithm and allow for the choice of TTM algorithm based on the relative sizes of dimensions as described in Sec. 6.3. To simplify the analysis, we provide an upper bound on the communication costs by assuming the Reduce version, which sends more messages, is used for each mode.

The $n$th mode Gram performs ${R}_{n}^{\varolessthan}I_{n}^{2}{I}_{n}^{\varogreaterthan}/P^{\oast}$ flops (exploiting the symmetry of the output), the $n$th mode eigenvalue computation requires $(10/3)I_{n}^{3}$ flops, and the $n$th mode TTM performs $2{R}_{n}^{\varolessthan}R_{n}I_{n}{I}_{n}^{\varogreaterthan}/P^{\oast}$ flops. Thus, the leading order terms in the flop costs are

[TABLE]

The communication cost of the $n$th mode (new) Gram computation is $\beta\cdot({R}_{n}^{\varolessthan}I_{n}{I}_{n}^{\varogreaterthan}/P^{\oast}\cdot((P_{n}{-}1)/P_{n})+I_{n}^{2})+\alpha\cdot(P_{n}{-}1+2\log P^{\oast})$. The communication cost of the $n$th mode (Reduce) TTM is $\beta\cdot({R}_{n}^{\varolessthan}R_{n}{I}_{n}^{\varogreaterthan}/{P}_{n}^{\oslash}\cdot((P_{n}{-}1)/P_{n}))+\alpha\cdot(2P_{n}\log P_{n})$. Thus, the leading order terms in the bandwidth costs, assuming we use the new Gram algorithm, are

[TABLE]

The leading order terms in the latency costs, conservatively assuming we use the Reduce TTM algorithm at each step, are

[TABLE]

The temporary memory required for the $n$th Gram computation is twice the size of the current local tensor data, which is ${R}_{n}^{\varolessthan}I_{n}{I}_{n}^{\varogreaterthan}/P^{\oast}$. Because this memory can be re-used and $R_{n}\leq I_{n}$ for each $n$, this cost is dominated by the first mode, requiring $2I^{\oast}/P^{\oast}$ words. The Gram computation also requires space for storing the output Gram matrix, which is of size $I_{n}^{2}/2$ and can be re-used across modes. The eigenvalue computation requires as much as $I_{n}^{2}$ extra memory if all eigenvectors are computed. The temporary memory required for TTM is guaranteed to be smaller than the input tensor, which is always bigger than the output tensor in the case of ST-HOSVD, so the temporary memory required for the $n$th TTM never exceeds the size required of the $n$th Gram computation. Thus, the leading order terms of the total temporary memory required on each processor is

[TABLE]

8. Optimized reconstruction

After the data has been compressed using ST-HOSVD, the user may wish to move the compressed data to another machine and reconstruct an approximation of the original data there. Since the full reconstructed data set would take up just as much space as the original, we provide the user an option to reconstruct sub-tensors of the original tensor. Thanks to the structure of the Tucker model, linear operations can be applied cheaply to each mode. For example, the user may want to downsample one of the spatial dimensions (mode 2) and select a single variable (mode 3) at a single timestep (mode 4) for the purposes of visualization. This can be accomplished through the following series of TTMs.

[TABLE]

where

[TABLE]

The five TTMs can be done in any order to obtain the result, but the Multi-TTM ordering can have a large effect on computational cost, memory footprint, and communication cost (in the parallel case). We demonstrate the effects of reconstruction mode ordering on run time and memory in Sec. 9.5. In the code, the user can specify the mode ordering or allow the software to automatically select the ordering that minimizes either computational cost or memory footprint. In the case of parallel reconstruction with linear operations applied to the factor matrices, each processor redundantly stores each ${\mathbf{C}}_{n}$ and computes ${{\mathbf{C}}_{n}}^{T}{\mathbf{U}}_{n}$ locally.

To determine the optimal mode ordering, the TuckerMPI code exhaustively searches over all $N!$ permutations to obtain the optimal ordering, but the optimal ordering can be determined in $O(N\log N)$ time by sorting with a specific comparator (Chakravarthy17, ). For example, consider the $N=2$ case, which corresponds to a product of 3 matrices, with input core matrix dimensions $R_{0}\times R_{1}$ and output subtensor dimensions $K_{0}\times K_{1}$. Computing the mode-0 product followed by the mode-1 product requires $K_{0}R_{0}R_{1}+K_{0}R_{1}K_{1}$ scalar multiplications, while computing the products in the opposite order requires $R_{0}R_{1}K_{1}+R_{0}K_{0}K_{1}$ scalar multiplications, so to minimize flops we order the products based on the comparison of these two costs.

More generally, let the input core tensor ${\bm{{\mathscr{{G}}}}}$ have dimensions ${R}_{0}^{\phantom{{}^{\prime}}}\times{R}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{R}_{N-1}^{\phantom{{}^{\prime}}}$ and the output subtensor ${\bm{{\mathscr{{Z}}}}}$ have dimensions ${K}_{0}^{\phantom{{}^{\prime}}}\times{K}_{1}^{\phantom{{}^{\prime}}}\times\cdots\times{K}_{N-1}^{\phantom{{}^{\prime}}}$. The key insight is that optimal ordering of all $N$ modes is the one such that every pair of modes is ordered correctly according to the $N=2$ case. That is, to minimize flops, mode $i$ should precede mode $j$ in the ordering if $K_{i}R_{i}R_{j}+K_{i}R_{j}K_{j}<R_{i}R_{j}K_{j}+R_{i}K_{i}K_{j}$. This can be shown, as argued by Chakravarthy (Chakravarthy17, ), by recognizing that this comparator yields a total ordering on modes and arguing by contradiction. Suppose the optimal ordering is not sorted by this comparator, then there exists two consecutive modes $n$ and $n+1$ that are out of order. Swapping the two consecutive modes will reduce the overall cost because the costs of the first $n-1$ TTMs and the last $N-n-1$ TTMs are equivalent, but the cost of the $n$ and $n+1$ TTMs are reduced by the property of the comparator, and we have a contradiction.

Similar arguments can be made for communication (bandwidth) cost and memory footprint. Given the bandwidth cost of Alg. 5, to minimize words moved, mode $i$ should precede mode $j$ if $R_{j}K_{i}(P_{i}-1)+K_{i}K_{j}(P_{j}-1)<R_{i}K_{j}(P_{j}-1)+K_{i}K_{j}(P_{i}-1)$. To minimize temporary memory, mode $i$ should precede mode $j$ if $R_{i}/K_{i}<R_{j}/K_{j}$.

9. Experimental Results

9.1. Experimental Platform

We run all experiments on Skybridge, a Sandia supercomputer consisting of 1,848 dual-socket 8-core Intel Sandy Bridge (2.6 GHz) compute nodes. Each node has 64 GB of memory, a peak flop rate of 332.8 GFLOPS (i.e., 20.8 GFLOPS per core), and the nodes are connected by an Infiniband interconnect. We use Intel compilers and the MKL for BLAS and LAPACK subroutines. We execute 16 MPI processes per node with 1 thread per process unless otherwise stated. Data files are stored on a Lustre file system, and Skybridge’s I/O is shared with other clusters. All reported timings in this section are averages over multiple runs. All reported memory requirements are computed by the application, using a wrapper around memory allocation calls, and does not include memory used by the operating system or by the MPI implementation.

9.2. Data Description

We consider two large-scale simulation datasets which were produced by S3D (CC+09, ), a massively parallel direct numerical simulation of compressible reacting flows, developed at Sandia National Laboratories. The datasets are as follows:

• •

SP: This 5-way data tensor is of size $500\times 500\times 500\times 11\times 400$ and corresponds to a cubic $500\times 500\times 500$ spatial grid for 11 variables over 400 time steps. Each time step requires 11 GB, so the entire dataset is 4.4 TB. The SP dataset is from the simulation of a 3D statistically steady planar turbulent premixed flame of methane-air combustion (KZCS16, ). The first 50 time steps (a 550 GB dataset) was used in previous work (AuBaKo16, ).

• •

JICF: This 5-way data tensor of size $1500\times 2080\times 1500\times 18\times 10$ comes from a $1500\times 2080\times 1500$ spatial grid with 18 variables over 10 time steps. Each time step requires 674 GB storage, so the entire dataset is 6.7 TB. The JICF dataset is from a jet in crossflow simulation, which is a canonical configuration for many combustion systems (LW+15, ).

The general experimental setup is described in Tab. 4. For each dataset, we use the same number of nodes for all experiments: 250 for SP and 350 JICF. Since each node has 64 GB of RAM and runs 16 threads/processes, storage of the full tensor requires a little more than $1/4$ of the memory per node (our data is stored in double precision). For each dataset, we consider three different processor grid configurations (A/B/C), as specified in Tab. 4b. They vary primarily in how the processors are distributed in the first three modes with scenario A having fewer processors in mode 0, scenario B being more evenly divided on modes 0–2, and scenario C having more processors on mode 0.

For testing and debugging, TuckerMPI also provides both sequential and parallel synthetic tensor generators for the users’ convenience. The user specifies the desired size and rank of a tensor, denoted by $(I_{0},I_{1},\dots,I_{N-1})$ and $(R_{0},R_{1},\dots,R_{N-1})$ respectively, and the amount of relative noise to be added, denoted $\eta$. We construct a core tensor ${\bm{{\mathscr{{G}}}}}$ with dimensions defined by $(R_{0},R_{1},\dots,R_{N-1})$; the numbers are drawn from a standard normal distribution. The factor matrices ${\mathbf{U}}_{0}\cdots{\mathbf{U}}_{N-1}$ are generated in a similar fashion, where ${\mathbf{U}}_{n}$ is $I_{n}\times R_{n}$. By performing the tensor times matrix mulitplications ${\bm{{\mathscr{{G}}}}}\times_{0}{\mathbf{U}}_{0}\times_{1}{\mathbf{U}}_{1}\cdots\times_{N-1}{\mathbf{U}}_{N-1}$, we obtain a tensor of size $(I_{0},I_{1},\dots,I_{N-1})$ with rank $(R_{0},R_{1},\dots,R_{N-1})$, which we refer to as ${\bm{{\mathscr{{M}}}}}$.

After obtaining the tensor ${\bm{{\mathscr{{M}}}}}$, we wish to add noise to it so that it is not exactly rank ${\bm{\mathbf{{R}}}}$. Let ${\bm{{\mathscr{{X}}}}}={\bm{{\mathscr{{M}}}}}+\eta\frac{\left\|{\bm{{\mathscr{{M}}}}}\right\|}{\left\|{\bm{{\mathscr{{N}}}}}\right\|}{\bm{{\mathscr{{N}}}}}$, where ${\bm{{\mathscr{{N}}}}}$ is a randomly generated tensor of noise whose values are also obtained from a standard normal distribution. This gives us $\frac{\left\|{\bm{{\mathscr{{X}}}}}-{\bm{{\mathscr{{M}}}}}\right\|}{\left\|{\bm{{\mathscr{{X}}}}}\right\|}\approx\eta$. Note that we do not explicitly construct and store ${\bm{{\mathscr{{N}}}}}$ due to its large size, so we approximate its norm using its expected value: $\left\|{\bm{{\mathscr{{N}}}}}\right\|\approx\sqrt{I^{\oast}}$.

9.3. Comparison of Gram Algorithms

We compare the two versions of the Gram algorithms described in Sec. 6.2: the old round robin variant from (AuBaKo16, ) and the new redistribution variant from Alg. 4. We use the two datasets and pick two representative modes, with corresponding experimental conditions detailed in Tab. 4. The results are shown in Fig. 7a. The new algorithm is up to 48 times faster than the old algorithm and never slower. The least speedup occurs when there is only one processor in the Gram mode fiber and no communication is performed in either case. The most speedup is when there are 100 processors in each fiber of the Gram mode. Ignoring the cost of the All-Reduce, the communication cost ratio between the old and new algorithms for a Gram operation in mode $n$ is $P_{n}$, which is an upper bound on the possible speedup. We see speedups of around 20-50% of that upper bound because of time that both algorithms spend on computation and the increased cost of the All-Reduce for the new algorithm. With respect to processor configurations, there is much less variance in the new algorithm than in the old. The timing for the new algorithm varies by no more than $3\times$, whereas the timing for the old algorithm varies between 12–62$\times$ depending on the processor configuration.

A breakdown of the run time for each of the experiments is shown in Fig. 7. We see that the main speedup comes from a reduction in communication of the tensor, but we also see a reduction in computation time. This is the result of the new algorithm making a single BLAS call rather than the old algorithm’s $P_{n}$ BLAS calls on smaller subproblems. The new algorithm does include some overhead for packing and unpacking the data, but it is negligible compared to the benefits of reduced communication.

Fig. 7a also reports the per-process memory requirement of each algorithm. Because the new algorithm requires re-packing the data and performing an All-to-All collective, it requires space for 2 extra copies of the local tensor data. The old algorithm requires temporary space for only one extra copy of the local tensor data, to perform its round robin exchange of data. The ratio of three total copies to two total copies yields the memory footprint ratio of 1.5 as reported in the table. When the number of processors in the mode of the Gram computation is one, no communication is necessary and no extra memory is required, so the memory footprint ratio is 1 in those cases.

9.4. ST-HOSVD

We analyze the parallel ST-HOSVD (Alg. 1) using the new redistribution version of the Gram algorithm. First, we consider its performance on two real-world datasets, varying the processor grid but not the number of processors. Second, we do strong and weak scaling studies varying the number of processors.

9.4.1. Compression of combustion data

We use TuckerMPI to compress the SP and JICF datasets described in Tab. 4a using the number of processors and grids specified in Tab. 4b. In terms of the pre-processing described in Sec. 3.5, there was no preprocessing on the SP data since it was already scaled, and max scaling was applied to the JICF data. For each dataset, we consider two relative error tolerances: 1e-2 and 1e-4, which we refer to as “High” and “Low” compression, respectively. The scenarios are summarized in Fig. 8a along with the compressed core size, total storage for the core and factor matrices, and overall compression ratio. For real-world datasets such as these, the compression potential depends on the amount of redundancy that is inherent in the data. For time-evolving simulations, there may be spatial regions with minimal change, and so these parts can be highly compressed. In our code, the user specifies the desired relative error tolerance ($\epsilon$), from which the level of compression is determined on the fly. (Alternatively, the software allows the user to specify the desired final core size, from which the final relative error is determined.) For our two datasets, the high-compression scenario yields reductions in size of 4–5 orders of magnitude. We note that our subject-matter experts have deemed the high compression datasets to be faithful representations that are scientifically useful. For instance, Fig. 9 shows a visual comparison between the original and compressed versions. Specifically, Fig. 9 shows a temperature isosurface at the 201st (middle) time point, used to track the boundaries of a flame during the simulation. The errors in the reconstructions of the compressed versions are unobservable in this visualization. As we see from Fig. 8, the resulting compressed data sets are small enough to be easily shared across high-speed networks and/or analyzed on workstations as described in Sec. 9.5.

Fig. 8 shows the run time results of ST-HOSVD using the different compression scenarios in Fig. 8a and processor grid configurations in Tab. 4b. Note that the degree of compression has no dependence whatsoever on the processors grid configuration; the configuration impacts only the run time. In the figure, the run times are broken down into color-coded segments corresponding to the three main kernels: Gram, Evecs, and TTM. As these are five-way tensors, there are five calls to each kernel. The times corresponding to mode zero are at the bottom and the most expensive. Since the tensor is reduced at each iteration (via the call to TTM), the calls become significantly cheaper so that the breakdown is hardly apparent past the first 2-3 iterations/modes. Comparing the “high” and “low” scenarios, the initial Gram computations are roughly equivalent, but all other computations in the 1e-2 case finish more quickly than the 1e-4 case because the data is compressed more drastically at each step. Comparing SP and JICF datasets, the Evecs computation is negligible for SP but is noticeable for JICF because its largest mode sizes are 3–4 times larger than for SP (the eigenvector computation is not parallelized). Comparing across processor grid configurations (A/B/C), loading more processors onto later modes generally improves run time, as the heavy communication steps are performed on data that can be orders of magnitude smaller than the initial tensor. The fastest of the processor grids are 3–6$\times$ faster than the slowest of the processor grids. This is not included in the figure, but the average time taken for preprocessing the JICF data (scaling each mode-3 slice by the inverse of the maximum entry) is 8 seconds and varied minimally across processor grids. This amounts to at most 15% of the run time of ST-HOSVD after the tensor is loaded in memory. The SP data in our experiments has already been scaled, so TuckerMPI does no pre-processing.

We note that this experiment does not consider alternative mode orderings, which can have a significant effect on run time. In the case of ST-HOSVD with a specified tolerance, the core tensor size is not known a priori, so it is not possible to pick an optimal ordering as discussed in Sec. 8. However, one can use some heuristics to pick a mode ordering. For example, starting with a mode whose processor grid dimension is 1 avoids communication of the tensor in the first Gram (which is typically the most expensive operation). Indeed, in this experiment, the fastest processor grid (A) has one processor in mode 0 for both SP and JICF data sets, but using the natural mode ordering 01234 is not necessarily optimal, even for that processor grid.

Fig. 8d shows the memory footprints and I/O times of the different cases. The input SP tensor requires about 1.1 GB per core on 4000 cores, and TuckerMPI requires about 3$\times$ that much memory to complete its computations, which is dominated by the memory required by the initial Gram computation. When a processor grid with only one processor in mode 0 is used, the memory footprint is only about 10-30% more than the initial data, depending on how much compression is achieved through the algorithm. The JICF data on 5600 cores requires about 1.2 GB per core, and again we see a memory footprint of about $3\times$ that amount. In terms of I/O time, the TuckerMPI code uses the MPI I/O interface for reading and writing binary files of multidimensional arrays. From the input results, we see a file reading bandwidth of 3–10 GB/sec., requiring $O(10)$ minutes to read the input tensors from disk. In comparison, computing the ST-HOSVD takes $O(1)$ minute per Fig. 8d. For instance, the SP-High-A scenario is 100$\times$ faster than reading the data from from disk. Writing to disk is much more expensive than reading. As a result, writing the compressed representation takes a disproportionate amount of time. We believe this is an artifact of the underlying MPI I/O implementation with the Lustre filesystem on Skybridge. One experiment, labeled ‘DNC’, did not complete due to an error within the MPI I/O implementation. Improving I/O performance for the low compression scenario is a topic of future work.

9.4.2. Strong scaling with synthetic data

In this section, we demonstrate the scalability of our code on a synthetic 4D tensor with dimensions $256\times 256\times 256\times 256$. The experimental setup is as follows. We fix the core dimensions to $32\times 32\times 32\times 32$ (4096$\times$ compression) so that all mode orderings are equivalent, and we scale from 1 node up to 128 nodes (2048 cores), using 1 MPI process per core. Since all the dimensions are the same, the mode ordering is irrelevant. We used the following processor grids: $1{\times}1{\times}2{\times}8$, $1{\times}1{\times}2{\times}16$, $1{\times}1{\times}4{\times}16$, $1{\times}1{\times}8{\times}16$, $1{\times}1{\times}16{\times}16$, $1{\times}2{\times}16{\times}16$, $1{\times}4{\times}16{\times}16$, $1{\times}4{\times}16{\times}32$, chosen heuristically to have fewer processors in the first and second modes to minimize communication in the early Gram computations, which dominate the run time.

The run times are reported in Fig. 10a. The input tensor is 32 GB in size, about half the memory available on a single node. The performance scales well up to 16 nodes (256 cores), achieving $9\times$ speedup over the single node, because the run time is dominated by the local Gram computation in the first node, which is perfectly parallelized. The degradation of performance after 16 nodes is caused in large part by the local computation in the first Gram computation (a single call to syrk) failing to scale perfectly. This is due to the local dimensions becoming too small and performance variability across nodes. We note that performance reported in (AuBaKo16, ) shows strong scaling of nearly the same algorithm on nearly the same problem to 256 nodes, but it was benchmarked on a different parallel computer.

9.4.3. Weak scaling with synthetic data

In this section, we demonstrate the weak scalability of our code on synthetic 4D tensors whose sizes scale with the number of processors so that the portion per processor remains constant. We ran our code on $2^{k}$ nodes (for $1\leq k\leq 5$), using 16 MPI processes per node, with a $1\times 1\times 4k^{2}\times 4k^{2}$ processor grid. We take a tensor of order $200k\times 200k\times 200k\times 200k$ and reduce it to a core size of $20k\times 20k\times 20k\times 20k$. This means that the local data on each node is approximately 0.8 GB. The processor grid was chosen as the best performing among $2k\times 2k\times 2k\times 2k$, $k\times k\times 4k\times 4k$ and $1\times k\times 4k\times 4k^{2}$. The largest speedup of the $1\times 1\times 4k^{2}\times 4k^{2}$ grid over the other grids was about 70%, and the $1\times k\times 4k\times 4k^{2}$ grid performed nearly as well.

Fig. 10b reports the performance in terms of GFLOPS per core. We observe that weak-scaling performance is generally preserved up to 625 nodes, with the code achieving 40-50% of peak performance throughout. In this experiment, the amount of local computation is held fixed, while the communication is increased with the number of processors. The weak scaling is possible because the run time remains dominated by computation. On 1 node ($k=1$), the tensor size is 12.8 GB, and on 625 nodes ($k=5$) it is 8 TB. Compared to the performance reported in (AuBaKo16, ), we observe similar overall performance relative to the architecture.

9.5. Reconstruction

Compressing data makes it cheaper to store and to transmit, but ultimately it needs to be reconstructed to be useful. We expect that most users will do only partial reconstructions. These can be used to visualize a portion of the data, extract summary statistics, etc. Moreover, this can oftentimes be done on a workstation rather than a parallel computer, making analysis much simpler. We can also reconstruct the entire dataset, which is primarily useful for comparison to the original dataset for quality control. In either case, reconstruction is a Multi-TTM operation as described in Sec. 8.

9.5.1. Partial Reconstruction

Our aim in this section is to show that we can do partial reconstructions on a laptop or workstation, so we run all experiments sequentially using only one MPI process (and thus only one node and one thread). We note that each node on Skybridge has 64 GB, and that the sequential execution has access to all 64 GB. We consider the SP dataset, as described in Sec. 9.2.

Visualization of Single Time Step

One partial reconstruction for visualization is to extract the entire $500\times 500\times 500$ grid corresponding to a single variable (out of 11) at a single timestep (out of 400). This results in a tensor of size $500\times 500\times 500\times 1\times 1$ which requires 1 GB of storage. This is the first step used to generate the visualizations of the compressed data in Fig. 9, for example.

As discussed in Sec. 8, the mode ordering is important in the reconstruction. We never want the intermediate size to be bigger than the larger of the input and output tensor sizes, but this can happen for the wrong mode ordering. Consider the high compression scenario and the intermediate tensor sizes that result as shown in Tab. 5. Using the ordering 01234 results in a 66 GB tensor after the third TTM in the reconstruction. But ordering 43120 never has an intermediate result larger than the final reconstruction. We stress that any ordering produces the same result (up to floating point error) and the difference is in the size of the intermediate results and the run time.

In Tab. 6, we report the maximum memory usage, the compute time (Multi-TTM), and the I/O times for different mode orders. The sizes of the core tensors on disk are 21 MB and 11 GB per Fig. 8a, and the size of the reconstructed output is 1 GB. The I/O time shows a bandwidth rate of about 1 GB/sec. In the high compression case, the partial reconstruction is larger than then the input; in contrast, the reverse is true for the low compression scenario. For both the high and low compression scenarios, the unique minimizer of the flops is mode order 43120, and it is one of the memory minimizers. Compared with the other orders benchmarked in the experiment, the optimal order runs as much as an order of magnitude faster, and some orders, including the straightforward 01234 order, fail due to out-of-memory errors.

Computing a Summary Statistic

As an example summary statistic for the SP data, we can compute the average fraction of carbon dioxide across all space and time in a few seconds. Mathematically, in the notation of Sec. 8, we set ${\mathbf{C}}_{0}$, ${\mathbf{C}}_{1}$, ${\mathbf{C}}_{2}$, and ${\mathbf{C}}_{4}$ to be all-ones vectors and ${\mathbf{C}}_{3}$ (the mode corresponding to variables) to be all zeros except for a one in the index that corresponds to carbon dioxide. Tab. 7 shows the results of this computation on a single node for both compression scenarios. The average fraction of carbon dioxide across all space and time is 0.0534 (calculated using the original data).111In preprocessing, the data was rescaled to have a maximum value of 0.5, and the mean based on the rescaled data was 0.3014. In comparison, the relative error of the averages computed from the compressed versions are 6.7e-6 (high compression) and 5.7e-8 (low compression). The memory required is only that of storing the compressed data, and the time is dominated by the cost of reading the core tensor from disk. For the high compression case, the computation took less than 1 second in total; for the low compression case, it required less than 16 seconds. Recall that working with the original data set requires a parallel computer and 100s of nodes just to read the data.

9.5.2. Full reconstruction of the SP data

As mentioned above, we do not expect users to employ full reconstruction very often, but it is useful as a diagnostic tool to check the quality of the approximation. Hence, we reconstructed the full SP dataset for both the high and low compression scenarios using a variety of mode orderings and processor configuration A ($1\times 1\times 40\times 1\times 100$) since it resulted in the fastest compression time. The results are reported in Tab. 8.

The I/O timings are averaged over the runs with different mode orderings. As the output is about 400 times larger than the larger of the two inputs, the overall time for the experiment is dominated by writing the output to disk, which takes almost an hour. (This supports the idea of avoiding reading and writing the full data sets to disk.) As compared with the bandwidth rate of one MPI process (see Tab. 6), the parallel I/O bandwidth rate is slower for the smaller inputs and about 50% faster for the larger output.

As discussed in Sec. 8, the computational and communication costs for the Multi-TTM, as well as the temporary memory footprint, are all dependent on the mode ordering of the individual TTMs. Furthermore, the mode ordering that minimizes computation need not be the same as the one that minimizes communication or memory. For the high compression scenario, mode order 41203 minimizes computation, 42103 minimizes communication, and 34120 minimizes memory. We also experiment with ordering 12034, which requires about 60% more temporary memory than the optimal orders and yields an out-of-memory error in the high compression scenario. The 34120 order yields minimum max memory usage for both scenarios, and 42103 is the fastest for both scenarios.

10. Conclusion

The Tucker tensor decomposition is useful for compression of many large-scale datasets because it uncovers latent low-dimensional structure. For multi-terabyte datasets that arise in direct numerical simulation of combustion reactions, we demonstrate that the Tucker decomposition can obtain five orders of magnitude in compression by reducing a 4.4 TB dataset to 21.5 MB. The compression time is an order of magnitude faster than simply reading the data. Austin, Bader, and Kolda (AuBaKo16, ) proposed the first parallel Tucker implementation. In this paper, we build upon their work, explaining the details of the parallel and serial algorithms as well as the local and distributed data layouts. Additionally, we have improved upon the algorithm in (AuBaKo16, ) and so can now run on much larger datasets. We also explain in detail how to reconstruct portions of the data, one of the most important benefits of the Tucker decomposition, without requiring any parallel resources. We show that it is possible to reconstruct portions of the data on a single processor in only a few seconds and using no more memory than the size of the input or output (whichever is larger).

Other parallel algorithms and implementations for computing Tucker decompositions of dense tensors using Higher-Order Orthogonal Iteration include (CC+17, ), which uses TTM-trees and performs dynamic regridding to optimize computation and communication costs, and (MS18, ), which uses TTM-trees and approximate information to update factor matrices in each iteration. Choi, Lui, and Chakaravarthy (CLC18, ) implement ST-HOSVD for dense tensors on a cluster equipped with GPUs, and they also use randomized SVD algorithms to improve run time. In the sparse case, there have also been several algorithmic innovations to improve the parallel performance of Higher-Order Orthogonal Iteration (KU16, ; SK17, ; OPLK18, ); these optimizations exploit the sparsity of the data tensor and avoid memory overhead of temporary dense data structures.

In future work, we hope to create a library for in situ compression that can compress at each time step in a simulation. Ideally, the simulation would never write the full dataset to disk but only compressed versions. Although these compressed versions cannot be used for restart in the event of a failure, they are a sort of “thumbnail” of the full simulation. This would save both time (for I/O) and disk space, not to mention making the sharing of data much easier.

We have not compared our approach to other methods of compression but hope to do so in future studies. Ballester-Ripoll, Lindstrom, and Pajarola (BaLiPa19, ) have recently compared Tucker compression (computed in serial) to the ZFP, SZ, and SQ compression methods and determined that this method “typically produces renderings that are already close to visually indistinguishable to the original data set.” Most other compression methods are focused on pointwise errors and achieve only about $O(10)$ compression and moreover divide up the space into blocks and compress the blocks individually. The Tucker approach we use here bounds the overall error and looks for large-scale patterns with the advantage being that it can get much higher compression, e.g., up to $O(10^{5})$. A comparison study would need to determine how to compare across different types of compression metrics (e.g., overall versus pointwise) and have a common large-scale dataset for the comparisons. We note that we can also potentially combine compression methods, using other methods to compress the core that results from the Tucker compression. Indeed, the “readme” file for the SZ Fast Error-Bounded Scientific Data Compressor lists Tucker compression as an option.222See https://github.com/disheng222/SZ/blob/master/README (ver. 12c9356), line 133.

Li et al. (LiGrPoCl15, ) evaluate the impact of wavelet compression (up to 512$\times$) for turbulent-flow data on visualization and analysis tasks, and such analysis would be interesting to consider also in the case of Tucker compression.

The data we compress in our study corresponds to a regular rectilinear grid. Another topic for future work is dealing with unstructured and/or non-rectilinear grids that are not neatly represented as a tensor. We also assume that the data is dense, but there are applications in data mining that have sparse tensors so this is another potential topic for study.

Acknowledgments

We would like to thank the anonymous referees for their comments and suggestions on the paper. We are indebted to Woody Austin for his work on the original prototype software on which this work is based. We thank Hemanth Kolla for creating Fig. 9 and his overall support of this work. We are also grateful to our colleagues Gavin Baker, Casey Battaglino, Cannada (Drew) Lewis, Jed Duersch, Alex Gorodetsky, and Prashant Rai for discussions about this software.

This material is based upon work supported in part by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA-0003525.