arXiv CSopen access

High-Performance Star-M SVD for Big Data Compression

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

cloud distributed-computing parallel-computing

distributed computing, parallel computing, cloud

High-Performance ⋆M SVD for Big Data Compression Md Taufique Hussain1 , Grey Ballard1 , Aditya Devarakonda1 , Srinivas Eswar2 , Naman Pesricha3 , and Vishwas Rao2

arXiv:2605.16058v1 [cs.DC] 15 May 2026

Wake Forest University, Winston-Salem, NC, USA 2 Argonne National Laboratory, Lemont, IL, USA 3 Indian Institute of Science, Bengaluru, India

Abstract In the era of big data, effectively compressing large datasets while performing complex mathematical operations is crucial. Tensor-based decomposition methods have shown superior compression capabilities with minimal loss of accuracy compared to traditional matrix methods. Under the ⋆M tensor framework, tensors can be decomposed in a matrix-mimetic way, including using the ⋆M SVD. This tensor SVD has optimality guarantees and has shown exceptional performance on specific types of data, but software implementations have been mostly limited to productivity-oriented languages. In this work, we present our development of a shared-memory parallel, high-performance solution designed to efficiently implement the underlying algorithms. This software will enable optimal compression of extensive scientific datasets, paving the way for enhanced data analysis and insights.

Keywords: parallel computing, high-performance computing, tensor decompositions

Introduction

Truncated tensor-based singular value decompositions based on the star-M (⋆M ) tensor-tensor product mimic the matrix SVD [26]. An input tensor can be decomposed into the product of three factor tensors with properties of orthogonality and structure that are analogous to the matrix case. Furthermore, this decomposition can be truncated for the purpose of data compression, for example, in a systematic way to determine the optimal representation given either an error tolerance or specified rank parameters. This is possible due to the existence of theoretical guarantees established for the ⋆M tensor algebra that generalize the Eckart-Young theorem [12] for low-rank matrix approximation. As described in more detail in Section 2, given 3-way input tensors, the ⋆M product is parametrized by an invertible matrix M and defined by a transformation of the mode-3 fibers of the inputs by M, independent matrix multiplications of the mode-3 slices, and a back transformation of the output mode-3 fibers by M−1 . The ⋆M truncated SVD (t-SVD) is thus based on slice-wise matrix SVDs after transformation by the matrix M, and it is particularly advantageous when the transformation yields many slices to be either low-rank or insignificant relative to other slices. The ⋆M t-SVD is not agnostic with respect to mode ordering, and the choice of M can have a substantial effect on the compressibility of a data set. Implementations of ⋆M algebra and the ⋆M t-SVD exist in MATLAB and Python—see [29, 31, 24], for example—but no high performance implementations have been developed. High performance, parallel software exists for many other tensor formats 1

including CANDECOMP/PARAFAC (CP) [17, 13, 37], Tucker [4, 22], and Tensor Train [2, 34]. The goal of this work is to design an efficient and scalable shared-memory parallel implementation of the fundamental algorithms for computing ⋆M t-SVD for general multiway data, enabling its use for datasets that were otherwise too large to handle using existing techniques. The core of the ⋆M t-SVD algorithm is a “batch” of independent matrix SVDs, all of the same dimension, which is ripe for parallelization. In fact, implementations of a batched SVD are emerging in standard libraries like MKL [19, Version 2024.1] and MAGMA [11, Version 2.10.0]. However, the (back-)transformations of the slices necessarily happens in the orthogonal direction, across slices, enabling the decoupling. We show in Section 3 that the transformations can be parallelized via efficient tensor-times-matrix operations that respect the data layout of the tensor in memory and exploits the interface to the batched BLAS [1]. Two primary algorithms exist for t-SVD: the first specifies a constant rank for the truncation of all slices, and the second specifies an error tolerance that guides the truncation to various ranks across slices. We show how to parallelize the first algorithm using the expert interface to LAPACK, and we explore computation-memory tradeoffs in the design of the second algorithm, which requires analysis of the computed slice-wise singular values to determine ranks. In Section 4, we consider data sets coming from the domains of climate reanalyses, computational fluid dynamics, and X-ray crystallography. We describe how mode ordering and choice of transformation matrix can affect both performance of the parallel implementations and the quality of the compressed representation. In the case of the climate data, we show the benefits of using tensor-based compression over the standard matrix-based approach used in that field. Our parallel design decisions enable efficient parallel scaling, and we show that on a platform with 64 cores, these datasets can be compressed up to 42×, 40×, and 30× faster than a single core, respectively. Our implementation (which will be made publicly available) enables users to observe speedups of up to 1000× over existing MATLAB code on a similar machine. Thus, this implementation allows domain scientists to analyze larger data sets in less time, enabling the ⋆M t-SVD to be applied as broadly as other tensor decomposition techniques.

Background

This section reviews the mathematical and algorithmic background relevant to our work. We first define notation for tensor and matrix-tensor operations, and define the ⋆M -product (Section 2.1). We then introduce the ⋆M t-SVD together with its full (t-SVDM-I) and adaptively truncated (tSVDM-II) variants, whose parallelization and evaluation is the main focus of our work (Section 2.2). Finally, we summarize the Eckart-Young-style optimality results that distinguish the ⋆M t-SVD from competing tensor decompositions (Section 2.3) and motivate the need for a shared-memory parallel implementation suited to large scientific datasets (Section 2.4).

t-SVDM-I

t-SVDM-II

rank ρi ≤ r varies per slice ρn discarded n

uniform rank r per slice n n

n ρn

ρ1

r Û ∈ Rm×r×n

p Ĝ ∈ Rr×p×n

ρ1 Û (variable rank)

p Ĝi (variable rank)

Figure 1: The output factor structure of t-SVDM-I (left) has uniform ranks across frontal slices which yield coarser error control, whereas t-SVDM-II (right) allows for variable frontal slice ranks ρi which can adapt to approximation requirements.

2.1

Notation

We use uppercase bold letters (e.g., M) to denote matrices and calligraphic letters (e.g., A) to denote tensors. A three-way tensor A ∈ Rm×p×n has n frontal slices A:,:,i ∈ Rm×p for i = 1, . . . , n. The mode-k unfolding A(k) arranges the mode-k fibers of A as columns of a matrix. The mode-k tensor times matrix product (TTM), written A ×k M for M ∈ RJ×nk , is defined by A ×k M (k) = M A(k) . The result has the same shape as A except that mode k has size J instead of nk . ⋆M -Product The ⋆M -product, introduced by Kilmer et al. [26], is a tensor-tensor multiplication parameterized by an invertible matrix M ∈ Rn×n . Given A ∈ Rm×p×n and B ∈ Rp×ℓ×n , the ⋆M -product C = A ⋆M B ∈ Rm×ℓ×n is defined in three steps. First, both operands are transformed to the M-domain by applying M along mode 3, Â = A ×3 M,

B̂ = B ×3 M.

Second, the n pairs of frontal slices are multiplied independently in the transform domain, Ĉ:,:,i = Â:,:,i B̂:,:,i ,

i = 1, . . . , n,

where each product is an ordinary matrix multiplication. Third, the result is mapped back to the original domain by applying M−1 along mode 3, C = Ĉ ×3 M−1 . When M is a nonzero scalar multiple of an orthogonal matrix, M−1 reduces to a (scaled) transpose, and the ⋆M -product inherits matrix-mimetic properties such as an associative and distributive algebra, a notion of orthogonality and identity tensors, and the Eckart-Young-style optimality results which are proved by Kilmer et al. [26]1 . 1

Extensions to 4th and higher order tensors are straightforward. Transformation to/from the M-domain occurs via TTMs with Mi /M−1 for i ≥ 3. A frontal slice is defined by fixing all indices except the first two for all the slice-wise i operations.

2.2

Tensor Singular Value Decomposition

In this section, we provide a brief overview of the mathematical background and algorithms that underpin our work. We focus on the tensor singular value decomposition (t-SVD) for compressing scientific datasets. The t-SVD can be defined in many ways, but we focus on the ⋆M t-SVD framework developed by Kilmer et al. [26], which is defined with respect to an invertible matrix M. Given an invertible matrix M ∈ Rn×n and a 3-way tensor A ∈ Rm×p×n , the ⋆M t-SVD of A is the factorization A = U ⋆M S ⋆M V, where U and V are orthogonal tensors and S is an f-diagonal tensor (i.e., diagonal frontal slices). The key computational step is transforming A to the M-domain by applying M to mode 3, Â = A ×3 M, after which the frontal slices Â:,:,i are mutually decoupled and each can be independently factored by a standard thin matrix SVD. The output factor tensors Û, Ŝ, V̂ are then mapped back to the original domain by applying M−1 to mode 3. The full ⋆M t-SVD is given in Algorithm 1, and the adaptive truncated variant with relative error tolerance ε is given in Algorithm 2 (see Figure 1). Choice of M. The transform matrix M is an important decision point in the ⋆M framework. Since M determines the basis used to decouple the frontal slices, it directly affects the singularvalue distribution and the resulting compression quality. Prior work has considered fixed analytic transforms (e.g. DCT and DFT), data-dependent transforms computed from A (e.g. singular vectors of the mode-3 unfolding), and learned transforms obtained by optimization. Kilmer et al. [26] survey several fixed, analytic choices that inherit the Eckart-Young optimality guarantee of the ⋆M t-SVD framework. The (normalized) discrete Fourier transform (DFT) recovers the original t-product of Kilmer and Martin and allows an O(n log n) application via the FFT, but yields complex-valued frontal slices in the transform domain. The discrete cosine transform (DCT) is a real-valued alternative that can also be applied with an O(n log n) algorithm. M can also be tailored to the input tensor. A simple data-driven choice discussed in [26] is to take M as the (transposed) matrix of left singular vectors of the mode-3 unfolding A(3) , which empirically yields better approximation at a given truncation rank than other fixed transforms. However, this data-driven approach comes at the cost of an upfront SVD on the unfolding. More recently, Newman and Keegan [31] have proposed a variable projection (VarPro) approach that jointly optimizes M and the truncated factors. The inner slicewise SVD problem is solved in closed form as a function of M and projected out, leaving a reduced nonlinear optimization problem. This produces transforms tailored to specific datasets and downstream tasks, at the cost of an iterative training step. Our work treats M as an arbitrary user-supplied invertible matrix, so any of the above choices can be utilized while leveraging the speedups achievable through shared-memory parallelization. Algorithm 1 Full t-SVDM of A Require: A ∈ Rm×p×n , invertible M ∈ Rn×n Ensure: U, S, V such that A = U ⋆M S ⋆M V⊺ 1: Â ← A ×3 M 2: for i = 1, . . . , n do 3: [Û:,:,i , Ŝ:,:,i , V̂:,:,i ] ← svd(Â:,:,i ) 4: end for

Algorithm 2 t-SVDM-II with relative error tolerance ε Require: Tensor A ∈ Rm×p×n , M ∈ Rn×n a nonzero multiple of an orthogonal matrix, relative error tolerance ε ∈ (0, 1) Ensure: Per-face truncation levels ρ = (ρ1 , . . . , ρn ) and truncated factors {Û:,1:ρi ,i , Ĝρi }ni=1 1: Compute the t-SVDM of A 2: v ← concatenate (Ŝj,j,i )2 for all i, j 3: v ← sort(v, ascending) Pk 4: Let w be the cumulative-sum vector, i.e., wk = ℓ=1 vℓ 2 2 5: Find the largest index J such that wJ /∥Ŝ∥F < ε 6: τ ← vJ ▷ largest singular value to be discarded 7: for i = 1, . . . , n do 8: ρi ← number of singular values of Â:,:,i greater than τ ⊺ 9: Keep only Û:,1:ρi ,i and form Ĝρi ← Ŝ1:ρi ,1:ρi ,i V̂:,1:ρi ,i 10: end for

2.3

Optimality of ⋆M t-SVD

The primary advantage of the ⋆M t-SVD is that it admits an Eckart-Young-like optimality guarantee whereas established tensor decompositions such as the CANDECOMP/PARAFAC (CP) [18, 8, 16], Tucker [36], Tensor-Train SVD [32], and the higher-order SVD (HOSVD) [9] families do not (see also the survey by Kolda and Bader [27]). For matrices, the Eckart-Young theorem states that truncating the SVD to its leading rank-k singular values and vectors yields the best rank-k approximation in the Frobenius (and 2-) norm [12, 30]. This single property is what makes the matrix SVD the primary tool for dimensionality reduction and lossy compression. Generalizing this guarantee to multiway data has been a long-standing open problem. For the CP decomposition, determining the tensor rank is NP-hard in general, and the factor matrices are not required to be orthogonal or even full rank [26]. As a result, no provably optimal truncation policy is known. Practical implementations often rely on iterative alternating-least-squares schemes that may converge slowly, are sensitive to initialization, and offer no global-optimality guarantees on the resulting approximation. The Tucker and HOSVD decompositions do produce factor matrices with orthonormal columns, but the core tensor is typically dense and does not yield an optimal low-rank approximation. Instead, the truncated HOSVD is only quasi-optimal : for a d-way tensor A, an HOSVD-truncated approximation BK satisfies √ ∥A − BK ∥F ≤ d ∥A − B⋆ ∥F , √ where B⋆ is the optimal multilinear-rank-k approximation [10]. The d factor compounds with tensor order, so the gap grows as the data become more high-dimensional. In contrast, when M ∈ Rn×n is a nonzero scalar multiple of an orthogonal matrix, Kilmer et al. [26] prove that both the t-SVDM-I (uniform-rank) and t-SVDM-II (error-tolerance, variable-rank) truncations are exactly Eckart-Young optimal in the Frobenius norm under the framework induced by ⋆M . This work showed that truncating the singular values of the transformed slices directly produces the best lower-rank representation. Combined with the fact that the slice-wise SVDs in the M-domain are decoupled (and easily parallelizable), the ⋆M framework provides an attractive option for compressing large scientific datasets where low-rank optimality and computational efficiency are important.

2.4

Shared-Memory Parallel ⋆M t-SVD

Despite the Eckart-Young-like optimality results, the practical impact of the ⋆M t-SVD is limited by the available software. The reference implementations accompanying [26] and the broader ⋆M work are written in productivity-oriented languages, such as MATLAB and pure Python. While these implementations are excellent for small scale experimentation with the ⋆M framework, they are not directly suitable for large-scale data compression tasks on modern multicore architectures. To our knowledge, no shared-memory parallel, high-performance implementation of the ⋆M t-SVD currently exists. This is in contrast to the software maturity of the matrix SVD (e.g. LAPACK and ScaLAPACK), CP, and Tucker decompositions, for which both shared-memory and distributedmemory parallel libraries are available [33, 35, 3, 5]. This work closes the parallel algorithms design and software gap by developing the pystarm software package, a shared-memory parallel, high-performance implementation of the ⋆M t-SVD that exposes parallelism in the tensor-times-matrix and the slicewise SVD operations. We utilize threaded BLAS and LAPACK kernels through Intel MKL, and provide a Python interface so that performance and productivity are simultaneously achievable for domain scientists.

Design and Implementation

We implement the two algorithms described in Section 2 in C++ using Intel’s Math Kernel Library (MKL) for dense linear algebra, with Python bindings provided via pybind11. Throughout, we restrict M to be orthonormal, so M−1 = M⊺ , which eliminates explicit matrix inversion. Recent work has extended the t-SVDM framework to other choices of M, including non-orthonormal transformations [25, 23]. We focus on orthonormal choices of M which provide optimality guarantees while reducing computational complexity when backward transformations of ⋆M t-SVD outputs are required. All tensor and matrix data are stored in column-major order. This choice is consequential: in column-major layout the m × p frontal slices of a 3-way tensor A ∈ Rm×p×n occupy contiguous, non-overlapping memory regions, so each slice can be addressed as a matrix without data movement or transposition. Both algorithms decompose into two main computational kernels—tensor-timesmatrix multiplication and slice-wise SVD—which we describe in turn before addressing the specifics of each variant. All microbenchmarks reported in this section for the two computational kernels were performed on a dual socket Intel Xeon CPU machine with 52 cores per socket.

3.1

Tensor-Times-Matrix

The TTM Â = A ×k M is the workhorse of both algorithms. For a 3-way tensor it applies M to mode 3 to enter the transform domain; for a d-way tensor it applies M to every mode k > 2 (with potentially a different transform per mode). The first two modes are never contracted by the algorithm, so we only need TTM kernels for modes 3, . . . , d. Among these, mode d is the last mode of the tensor and modes 3, . . . , d − 1 are middle modes; the structure of the computation differs in the two cases. For Q a TTM on the last mode of a d-way tensor, the column-major layout exposes the tensor as a ⊺ nd ×J produces the correctly contracted single i<d ni × nd matrix; right-multiplying it by M ∈ R result and is a single level-3 BLAS GEMM call [6]. For a 3-way tensor this is the only TTM that arises (k = d = 3), so MKL’s multi-threaded GEMM handles all parallelism automatically. For a TTM on a middle mode k with 2 < k < d, the mode-k fibers Q are not contiguous,Qso the operation must be expressed as a collection of GEMM calls. Let Mk = i<k ni and Pk = i>k ni .

The tensor can be viewed as Pk contiguous submatrices of shape Mk × nk , each to be rightmultiplied by M⊺ ∈ Rnk ×J to yield an Mk × J result [6]. Middle-mode TTMs only arise once the tensor has more than three modes, but then they dominate: in the 6-way ncep-air-6 dataset (Section 4.1), three of the four required TTMs (modes 3, 4, and 5) are middle-mode contractions. The efficient parallelization of these middle-mode TTMs is therefore an important design choice for a high-performance shared-memory implementation of the ⋆M t-SVD. 3.1.1

Multithreaded TTM Design Choices

We identify three primary choices to implement multithreaded TTM for the middle modes. Choice 1: Batched GEMM We can use MKL’s batched strided GEMM interface, in which a single MKL call performs all Pk multiplications such that: the input stride between successive submatrices is Mk nk , the output stride is Mk J, and the stride for M is zero so the same matrix is reused across all batches. This allows MKL’s scheduling layer to optimize across the batch, reducing per-call overhead compared to an equivalent loop over individual DGEMM calls. For the mode-3 TTM that appears in both algorithms on a 3-way tensor, Pk = 1 and Mk = mp, so the batched call degenerates to a single DGEMM that right-multiplies the (mp) × n matrix view of the tensor by M⊺ ∈ Rn×J , yielding an (mp) × J result. We call this the batched variant of TTM. Choice 2: Loop over multithreaded GEMM We can call the individual GEMM operations one after another sequentially and use optimized multithreaded GEMM kernels for each of the individual GEMMs. We call this the loop variant of TTM. Choice 3: Parallel for with sequential GEMM We can lift the parallelism from MKL multithreading to OpenMP multithreading by leveraging a parallel for over the batch of GEMM and use sequential MKL GEMM calls. These calls can be issued simultaneously without synchronization. We call this the parfor variant of TTM. 3.1.2

Performance of three TTM designs

We compare the performance of the three design choices of TTM on the ncep-air-6 tensor discussed in Section 4.1. It is a 6-way tensor and exercises both the middle-mode (modes 3, 4, 5) and last-mode (mode 6) cases that the algorithm actually invokes. We show the comparison in Figure 2. We observe that the batched variant works best independent of which mode the TTM operation is performed on. The parfor variant suffers for the last mode. For TTM on the last mode, the parfor variant reduces to a single sequential GEMM call, hence it does not scale at all with increasing thread count. The loop variant suffers for the first few modes. As explained above, TTM on a middle mode reduces to multiple GEMM calls. For the first few modes, the loop variant reduces to multiple calls of multithreaded GEMM each having a small workload. Hence, the scalability suffers at higher thread counts. For the remainder of this paper we report experiments with only the batched variant of TTM as it is the clear winner among all three design choices.

TTM Benchmarking

mode 3 (dimension=17)

time (s)

101

time (s)

ncep-air-6

100 1

threads mode 4 (dimension=4)

100

threads mode 5 (dimension=365)

101 100

threads mode 6 (dimension=10)

100

threads

batched

loop

parfor

Figure 2: TTM performance of three variants (batched vs. loop vs. parfor ) for the ncep-air-6 dataset across thread counts and TTM modes.

3.2

Slice-wise SVD

After transforming to the M-domain, both algorithms require computing the matrix SVD of each frontal slice Â:,:,i independently. Because LAPACK’s SVD routine overwrites its input, each slice is materialized as a fresh copy before factorization. We support two factorization modes. The first computes the full thin SVD of each slice, yielding factor tensors Û ∈ Rm×r×n and V̂ ∈ Rp×r×n (where r = min(m, p)) together with a matrix of singular values. The second computes only the leading k singular triplets via a truncated SVD, avoiding the cost of forming the full factorization when k ≪ r. 3.2.1

Multithreaded Design Choices

The n slice SVDs are independent, presenting two natural parallelization strategies. Choice 1: Sequential slices, parallel SVD Each LAPACK call is itself multi-threaded, exploiting intra-slice parallelism through MKL’s internal threading. Choice 2: Parallel slices, sequential SVD In the second, a parallel loop distributes slices across threads, with each thread calling a single-threaded LAPACK routine. The first strategy is preferable when the slices are large relative to the number of slices (m, p ≫ n), so each SVD has enough work to keep many threads busy; the second is preferable when there are many small slices (n ≫ m, p), where inter-slice parallelism offers more opportunity than intra-slice parallelism.

Slicewise SVD Benchmarking

time (s)

ncep-air-6 (248,200 slices of 73×144) 103 1.0x 102

102 1.0x

time (s)

4.0x

8.0x

101 1

101 1.0x

2.0x

4.0x

7.9x

31.2x

36.3x

15.7x

30.2x

46.3x

28.0x

26.8x

threads xray (400 slices of 300×400)

2.0x

3.9x

7.7x

100 1

16.0x

threads cfd (6,072 slices of 253×253)

101 100

time (s)

2.0x

threads

15.2x

parallel slices, sequential SVD sequential slices, parallel SVD Figure 3: Slice-wise SVD wall time (parallel slices, sequential SVD vs. sequential slices, parallel SVD) for the NCEP-Air 6-way, CFD, and X-ray Crystallography datasets. 3.2.2

Performance of the designs

In Figure 3 we show a microbenchmark demonstrating the performance of the two design choices of slicewise SVD on the three practical tensors described in Section 4. We observe that the parallelslices-with-sequential-svd variant gives almost linear scalability over increasing thread count. In contrast, the sequential-slices-with-parallel-svd variant does not scale. We believe that this is due to the small sizes of the SVD where parallel overhead dominates. Given the complexity of parallelizing the SVD computation, we expect the parallel-slices-with-sequential-svd to yield better performance for most applications.

3.3

Full t-SVDM and t-SVDM-I

Full t-SVDM Algorithm 1 shows the full t-SVDM algorithm. For simplicity, the algorithm is shown for a three-way tensor. This variant computes all singular values and vectors of each frontal slice. In contrast, t-SVDM-I computes a truncated SVD with fixed rank-r for each frontal slice. This variant enables predictable performance in each SVD call and reduces the overhead of computing variable rank SVDs and the complexity of determining ranks adaptively. The full SVD uses the dgesvd LAPACK call, whereas the truncated SVD calls dgesvdx.

3.4

t-SVDM-II

Tensor shape: [253, 253, 253, 4, 6] Time (s)

4 3 2 1 0

0.001

0.01

0.1

0.2

0.001

0.01

0.1

0.2

Stage 1 Threshold

Stage 2 truncate

Time (s)

Tolerance Tensor shape: [1500, 800, 20, 4, 6]

10 5 0

Tolerance

compute-efficient memory-efficient

Figure 4: Breakdown times for the different t-SVDM-II strategies. Depending on the time and memory constraints, any of the three methods can be preferred. The adaptive truncated variant (Algorithm 2) introduces per-slice truncation ranks ρi that vary across slices and are not known until after the singular values have been inspected globally. This data dependency rules out the simple composition used for t-SVDM-I and motivates three strategies with different trade-offs between computational redundancy and memory efficiency.

3.4.1

Multithreaded Design Choices

All strategies consist of three stages: (1) compute the singular values of all slices, (2) determine the threshold τ , and (3) compute the final factorization. Choice 1: Truncate The truncation strategy is a straightforward way to compute the t-SVDM-II factorization. We compute the full t-SVDM (Algorithm 1) factorization, determine the threshold τ from all n · r singular values, and discard the singular vectors corresponding to values below τ . This is algorithmically simple and reuses the same slice-wise SVD kernel without modification, but pays the full cost of computing all singular vectors only to discard most of them when compression is aggressive. The truncate variant consumes the most intermediate memory of all three strategies storing up to 3× the memory footprint of the input data. Choice 2: Compute-efficient. In the second strategy, a first pass computes only the singular values of each slice, without forming any singular vectors. The SVD algorithm consists of multiple steps (a possible QR factorization and bidiagonal reduction) before computing the singular values and vectors. We cache these intermediate entities. The global threshold τ and per-slice ranks ρi are determined from this pass. A second pass then computes only the leading ρi singular triplets for each slice, via a truncated SVD. We utilize the cached intermediaries to speed up the second SVD computation. All output buffers are formed to the exact required ranks before the second pass. The compressed factorization is stored in a variable-rank format in which the i-th slice’s factor matrices occupy exactly m × ρi and p × ρi entries with no padding, as opposed to the dense m × r × n and p × r × n tensors produced by t-SVDM-I. This representation accurately reflects the storage cost of the compressed result and avoids allocating memory for singular components that are not retained. When target ranks are small, the savings in the second stage outweigh the cost of reading each slice twice. This compute-efficient version of t-SVDM-II reduces the memory footprint by half when compared to the truncate strategy. It also performs the fewest floating-point operations of all three strategies for most cases. Choice 3: Memory-efficient Finally, we consider the most memory efficient strategy, where we do not save any intermediaries from the preceding strategy. The second SVD computation could possibly take longer than the compute-efficient strategy but we also save on the memory bandwidth needed in moving and copying the intermediaries. This strategy does not use any more space than that needed for holding the input and output buffers in memory. We call it the memory-efficient variant. 3.4.2

Performance of the three variants

Figure 4 shows the runtime breakdown of the different strategies on two data tensors for various relative-error tolerances. All times are averaged over 5 runs. The first tensor is the CFD data (253×253×253×4×6, ≈ 3 GiB) and the second is a larger tensor (1500×800×20×4×6, ≈ 4.3 GiB) with a similar singular value distribution. No clear winner exists among the three strategies. The truncate strategy is most insensitive to the tolerance parameter as it computes the full t-SVDM every time. Its runtime is fairly consistent over all tolerances with the bulk of the time in the Stage 1 SVD. The compute-efficient and memory-efficient strategies show more differing behavior. The Stage 2 times for compute-efficient are always lower than those of memory-efficient, but the gain may not always be sufficient to overcome the time spent in moving the cached intermediaries in memory during the first stage. The Stage 1 SVD calls the dgesvd kernel which employs the implicit 12

zero-shift QR algorithm for its bidiagonal SVD (dbdsqr), while the truncated SVD (dgesvdx) calls the bisection plus inverse iteration algorithm (dbdsvdx). Our studies showed an 8× performance difference between these methods, with dbdsqr outperforming dbdsvdx. This explains the relative competitiveness of the truncate strategy. Since we are concerned with data compression, we will use the memory-efficient strategy for all t-SVDM-II benchmarks.

Application

In this section we demonstrate the practical utility of our work on three scientific datasets that are prohibitively large for currently available t-SVDM-I and t-SVDM-II implementations: atmospheric reanalysis data (Section 4.1), computational fluid dynamics data (Section 4.2), and X-ray diffuse scattering data (Section 4.3). While tensor decompositions have many applications, we focus on compression as it provides a direct and quantifiable measure of the quality of the decomposition. For each dataset we evaluate compression quality and strong scaling performance.

4.1 4.1.1

Climate Data (NCEP-Air) Dataset

We use 10 years (1948–1957) of air temperature data from the NCEP/NCAR Reanalysis [21] at all 17 standard pressure levels, representing altitude. We represent the data as a 6-way tensor of shape 73 × 144 × 17 × 4 × 365 × 10 (ncep-air-6), where the first three modes are latitude, longitude, and pressure level, and the time axis is split into time-of-day (4), day-of-year (365), and year (10). February 29 observations are dropped from leap years to produce a uniform 365-day calendar, giving approximately 19.46 GB as float64. Compared to a flat 4-way representation (ncep-air-4, shape 73 × 144 × 17 × 14,612), in our experiments we observed that this reshaping achieves nearly identical compression quality while reducing TTM FLOPs by approximately 40×, as each time-mode GEMM operates on a much smaller dimension. 4.1.2

Compression Quality

Figure 5 compares compression ratio across methods on the ncep-air-6 tensor. t-SVDM-II-DCT achieves higher compression than t-SVDM-I-DCT at every error level, as per-slice variable truncation allows t-SVDM-II to discard many slices entirely. Both tensor methods outperform EOF [15], the standard matrix-based compression in climate science, consistent with the theoretical results of [26]. The DCT transform further outperforms the identity transform, as it concentrates energy into fewer singular values per slice. For the remaining experiments we focus on comparing t-SVDM-I and t-SVDM-II with the DCT transformation.

ncep-air-6: compression ratio

compression ratio

107

t-SVDM-I-DCT t-SVDM-II-DCT t-SVDM-II-Identity EOF

105 103 101

10 4

10 3

relative error

10 2

Figure 5: Compression ratio of the ncep-air-6 tensor for different algorithms. 4.1.3

Strong Scaling

Figure 6 shows strong scaling of t-SVDM-I (k = 1, compression ratio ≈ 48.2×) and t-SVDM-II (tol=0.005, compression ratio ≈ 48.9×) on the ncep-air-6 tensor. t-SVDM-I achieves 42× speedup from 1 to 64 threads (147s → 3.5s) and t-SVDM-II achieves 27× (160s → 5.9s). For t-SVDM-II, singular value computation dominates at low thread counts, while singular vector computation is relatively cheap — due to high compressibility, most slices are dropped entirely after thresholding.

number of threads

t-SVDM-I (k=1) t-SVDM-II (tol=0.005) transform TTM slicewise SVD

42.6× 27.3×

28.7× 20.4×

15.1× 12.5×

Strong Scaling 7.0× 6.4×

3.8× 3.6×

2.0× 1.9×

runtime (seconds)

ncep-air-6

slicewise SVD (values) thresholds slicewise SVD (vectors)

Figure 6: Strong scaling of both t-SVDM-I and t-SVDM-II for the ncep-air-6 data with DCT transformation. t-SVDM-II is run with tolerance 0.005 (compression ratio ≈48.9), and t-SVDM-I is run with rank k = 1 (compression ratio ≈48.2).

4.1.4

Extreme Event Preservation

EOF vs. t-SVDM-II at extremes. A known limitation of EOF-based compression in climate science is that it can smooth extreme weather events [20]. Figure 7 compares the median pointwise relative error at extreme temperature events (top and bottom 0.5% of the time series at each grid point, at 850 hPa — a pressure level that approximates surface air temperature) for EOF and t-SVDM-II-DCT at a matched compression ratio of approximately 320×. For each spatial grid point, the error is computed over all time steps in the dataset and aggregated as the median across extreme events. t-SVDM-II-DCT achieves lower reconstruction error at extremes across most of the globe, demonstrating that tensor-based compression better preserves physically significant rare events.

Reconstruction error at extreme events (850 hPa) [1948 1957] EOF

median relative error at extremes

0.045 0.040 0.035

0.025

tsvdmii

median relative error at extremes

Relative error

0.030

0.020 0.015 0.010 0.005

Figure 7: Median pointwise relative error at extreme temperature events (850 hPa) for EOF (top) and t-SVDM-II-DCT (bottom), at a matched compression ratio of ∼320×. Lower is better.

4.2 4.2.1

CFD Dataset

We use a computational fluid dynamics (CFD) dataset generated from a direct numerical simulation of the Taylor–Green Vortex (TGV) [7], a standard benchmark for turbulent flow, solved using NekRS [14]. The simulation data is assembled into a 5-way tensor of shape 253 × 253 × 253 × 4 × 6 (approximately 3 GiB as float64), where the first three modes index the structured spatial grid, 15

the fourth mode indexes four solution fields (velocity components u, v, w and pressure p), and the last mode indexes six flow snapshots. A visualization of one snapshot is shown in Figure 8.

Figure 8: Iso-surfaces of the z-component of vorticity, ωz = ∂v/∂x − ∂u/∂y, for the Taylor–Green Vortex flow: original (top-left), compressed reconstruction (bottom-left), and pointwise error (topright). Reconstruction is via t-SVDM-II-DCT at tolerance 10−1 (compression ratio ≈ 16×). Red and blue denote positive and negative values. 4.2.2

Compression Quality

Figure 9 shows the compression ratio of t-SVDM-I and t-SVDM-II on the cfd tensor. At tight tolerances (relative error ≤ 0.05), the two algorithms achieve comparable compression ratios — for example, at relative error ≈ 0.05, t-SVDM-I (k = 20) achieves 6.3× while t-SVDM-II achieves 5.4×. However, as the tolerance is relaxed, t-SVDM-II pulls significantly ahead: at relative error ≈ 0.1, t-SVDM-II achieves 16× compression compared to 6.3× for t-SVDM-I, and the gap grows further at higher tolerances. Figure 8 shows the reconstructed vorticity field at tolerance 10−1 , confirming that the compressed reconstruction closely matches the original.

compression ratio

cfd: compression ratio t-SVDM-I-DCT t-SVDM-II-DCT

103 102 101 100 10 4

10 3

10 1

10 2

relative error

Figure 9: Compression ratio of the cfd tensor for different algorithms. 4.2.3

Strong Scaling

Figure 10 shows strong scaling of t-SVDM-I (k = 20, compression ratio ≈ 6.3×) and t-SVDM-II (tol=0.05, compression ratio ≈ 5.4×) on the cfd tensor. t-SVDM-I achieves a 39× speedup from 1 to 64 threads (43s → 1.1s), and t-SVDM-II achieves 34× (61s → 1.8s). t-SVDM-II splits the SVD into two passes (singular values first, then singular vectors for retained singular values), which results in slightly weaker scaling. Unlike the ncep-air-6 tensor, singular vector computation accounts for a significant fraction of t-SVDM-II runtime — at 1 thread it takes ≈ 39s, comparable to the full SVD cost of t-SVDM-I.

39.5× 34.0×

27.5× 23.4×

2.0

14.6× 12.7×

8.0

7.6× 7.3×

32.0

Strong Scaling

3.8× 3.7×

2.0× 1.9×

runtime (seconds)

CFD

0.5 1

number of threads

t-SVDM-I (k=20) t-SVDM-II (tol=0.05) transform TTM slicewise SVD

slicewise SVD (values) thresholds slicewise SVD (vectors)

Figure 10: Strong scaling of both t-SVDM-I and t-SVDM-II for the cfd data with DCT transformation. t-SVDM-II is run with tolerance 0.05 (compression ratio ≈5.4), and t-SVDM-I is run with rank k = 20 (compression ratio ≈6.3).

4.3 4.3.1

X-ray Crystallography Dataset

We use a 3D X-ray diffuse scattering dataset measured on a MoVO2 single crystal [28], stored as a 3-way tensor of shape 300 × 400 × 400 in float64 (approximately 366 MB), where each mode corresponds to one of the three orthogonal axes of the volumetric reciprocal-space measurement. A representative slice through the volume is shown in Figure 11.

Figure 11: Two-dimensional slice through the X-ray diffuse-scattering volume: original (top-left), compressed reconstruction (bottom-left), and pointwise absolute error (top-right). Reconstruction is via t-SVDM-II-DCT at tolerance 10−2 (compression ratio ≈ 1.51×). A shared logarithmic colormap is used across all panels. 4.3.2

Compression Quality

Figure 12 shows the compression ratio of t-SVDM-I and t-SVDM-II on the xray tensor. t-SVDM-II achieves higher compression than t-SVDM-I at all error levels. Figure 11 shows the reconstructed slice at tolerance 10−2 , confirming that the reconstruction closely matches the original.

compression ratio

xray: compression ratio t-SVDM-I-DCT t-SVDM-II-DCT

102 101 100

10 4

10 3

10 1

10 2

relative error

Figure 12: Compression ratio of the xray tensor for different algorithms. 4.3.3

Strong Scaling

Figure 13 shows strong scaling of t-SVDM-I (k = 150, compression ratio ≈ 1.14×) and t-SVDM-II (tol=10−4 , compression ratio ≈ 1.19×) on the xray tensor. Both algorithms scale similarly — t-SVDM-I achieves 28× speedup from 1 to 64 threads (16.7s → 0.59s) and t-SVDM-II achieves 30× (18.7s → 0.63s). In both cases, slicewise SVD dominates the runtime.

23.6× 23.8×

28.2× 29.9×

1.0

14.2× 14.0×

7.5× 7.4×

4.0

Strong Scaling

3.8× 3.9×

16.0

1.9× 2.0×

runtime (seconds)

X-ray

0.2 1

number of threads

t-SVDM-I (k=150) t-SVDM-II (tol=0.0001) transform TTM slicewise SVD

slicewise SVD (values) thresholds slicewise SVD (vectors)

Figure 13: Strong scaling of both t-SVDM-I and t-SVDM-II for the xray data with DCT transformation. t-SVDM-II is run with tolerance 0.0001 (compression ratio ≈1.19), and t-SVDM-I is run with rank k = 150 (compression ratio ≈1.14).

Conclusions and Future Work

We present the first high-performance implementations of tensor decompositions in the ⋆M -framework. Our shared-memory parallel implementation takes advantage of the recent batch/strided interfaces developed for these linear algebra kernels. We achieve good scaling up to 64 threads on three real-world datasets from different scientific applications. Our case study on the climate dataset also highlights the utility of maintaining the multidimensionality of the dataset by preserving physically relevant rare events. Our initial investigations open up many future directions for research. Supporting GPUs in the library is an immediate next step. The embarrassingly parallel nature of the slice-wise operations would greatly benefit from this development. The irregular data access pattern due to varying ranks in t-SVDM-II also poses some difficulty. Developing distributed-memory parallel algorithms is a more involved next step. Overcoming the difference in data layout preferences for TTM and SVD stage without an explicit transpose operation is an immediate concern. Finally, we utilize explicit matrix multiplications in our TTM stage. Supporting operator transforms, for transforming the modes, is another avenue for future work.

References [1] Ahmad Abdelfattah, Timothy Costa, Jack Dongarra, Mark Gates, Azzam Haidar, Sven Hammarling, Nicholas J. Higham, Jakub Kurzak, Piotr Luszczek, Stanimire Tomov, and Mawussi Zounon. A set of batched basic linear algebra subprograms and lapack routines. ACM Trans. Math. Softw., 47(3), June 2021. [2] Hussam Al Daas, Grey Ballard, and Peter Benner. Parallel algorithms for tensor train arithmetic. SIAM Journal on Scientific Computing, 44(1):C25–C53, 2022. [3] Woody Austin, Grey Ballard, and Tamara G. Kolda. Parallel tensor compression for large-scale scientific data. In 2016 IEEE International Parallel and Distributed Processing Symposium (IPDPS), pages 912–922, 2016. [4] Grey Ballard, Alicia Klinvex, and Tamara G. Kolda. TuckerMPI: A parallel C++/MPI software package for large-scale data compression via the tucker tensor decomposition. ACM Transactions on Mathematical Software, 46(2), June 2020. [5] Grey Ballard, Alicia Klinvex, and Tamara G. Kolda. TuckerMPI: A parallel C++/MPI software package for large-scale data compression via the Tucker tensor decomposition. ACM Transactions on Mathematical Software, 46(2):1–31, 2020. [6] Grey Ballard and Tamara G. Kolda. Tensor Decompositions for Data Science. Cambridge University Press, 2025. [7] Shivam Barwey, Pinaki Pal, Saumil Patel, Riccardo Balin, Bethany Lusch, Venkatram Vishwanath, Romit Maulik, and Ramesh Balakrishnan. Mesh-based super-resolution of fluid flows with multiscale graph neural networks. Computer Methods in Applied Mechanics and Engineering, 443:118072, 2025. [8] J. Douglas Carroll and Jih-Jie Chang. Analysis of individual differences in multidimensional scaling via an N-way generalization of “Eckart-Young” decomposition. Psychometrika, 35(3):283– 319, 1970. 20

[9] Lieven De Lathauwer, Bart De Moor, and Joos Vandewalle. A multilinear singular value decomposition. SIAM Journal on Matrix Analysis and Applications, 21(4):1253–1278, 2000. [10] Lieven De Lathauwer, Bart De Moor, and Joos Vandewalle. On the best rank-1 and rank(r1 , r2 , . . . , rn ) approximation of higher-order tensors. SIAM Journal on Matrix Analysis and Applications, 21(4):1324–1342, 2000. [11] Jack Dongarra, Mark Gates, Azzam Haidar, Jakub Kurzak, Piotr Luszczek, Stanimire Tomov, and Ichitaro Yamazaki. Accelerating numerical dense linear algebra calculations with gpus. Numerical Computations with GPUs, pages 1–26, 2014. [12] Carl Eckart and Gale Young. The approximation of one matrix by another of lower rank. Psychometrika, 1(3):211–218, 1936. [13] Srinivas Eswar, Koby Hayashi, Grey Ballard, Ramakrishnan Kannan, Michael A. Matheson, and Haesun Park. PLANC: Parallel low-rank approximation with nonnegativity constraints. ACM Transactions on Mathematical Software, 47(3), June 2021. [14] Paul Fischer, Stefan Kerkemeier, Misun Min, Yu-Hsiang Lan, Malachi Phillips, Thilina Rathnayake, Elia Merzari, Ananias Tomboulides, Ali Karakus, Noel Chalmers, and Tim Warburton. NekRS, a GPU-accelerated spectral element Navier–Stokes solver. Parallel Computing, 114:102982, 2022. [15] A. Hannachi, I. T. Jolliffe, and D. B. Stephenson. Empirical Orthogonal Functions and Related Techniques in Atmospheric Science: A Review. International Journal of Climatology, 27(9):1119–1152, 2007. [16] Richard A. Harshman. Foundations of the PARAFAC procedure: Models and conditions for an “explanatory” multi-modal factor analysis. UCLA Working Papers in Phonetics, 16:1–84, 1970. [17] Ahmed E. Helal, Jan Laukemann, Fabio Checconi, Jesmin Jahan Tithi, Teresa Ranadive, Fabrizio Petrini, and Jeewhan Choi. ALTO: adaptive linearized storage of sparse tensors. In Proceedings of the 35th ACM International Conference on Supercomputing, ICS ’21, page 404?416, New York, NY, USA, 2021. Association for Computing Machinery. [18] Frank L. Hitchcock. The expression of a tensor or a polyadic as a sum of products. Journal of Mathematics and Physics, 6(1–4):164–189, 1927. [19] Intel Corporation. Developer Reference for Intel® oneAPI Math Kernel Library, 2026. [20] Yujing Jiang, Daniel Cooley, and Michael F. Wehner. Principal component analysis for extremes and application to U.S. precipitation. Journal of Climate, 33(15), 2020. [21] E. Kalnay, M. Kanamitsu, R. Kistler, W. Collins, D. Deaven, L. Gandin, M. Iredell, S. Saha, G. White, J. Woollen, Y. Zhu, M. Chelliah, W. Ebisuzaki, W. Higgins, J. Janowiak, K. C. Mo, C. Ropelewski, J. Wang, A. Leetmaa, R. Reynolds, R. Jenne, and D. Joseph. The NCEP/NCAR 40-year reanalysis project. Bulletin of the American Meteorological Society, 77(3):437–472, 1996. [22] O. Kaya and B. Uçar. High performance parallel algorithms for the Tucker decomposition of sparse tensors. In 45th International Conference on Parallel Processing (ICPP ’16), pages 103–112, 2016. 21

[23] Katherine Keegan and Elizabeth Newman. Projected tensor-tensor products for efficient computation of optimal multiway data representations. Linear Algebra and its Applications, 729:100–147, 2025. [24] Katherine Keegan, Tanvi Vishwanath, and Yihua Xu. A tensor svd-based classification algorithm applied to fmri data, 2021. [25] Eric Kernfeld, Misha Kilmer, and Shuchin Aeron. Tensor-tensor products with invertible linear transforms. Linear Algebra and its Applications, 485:545–570, 2015. [26] Misha E Kilmer, Lior Horesh, Haim Avron, and Elizabeth Newman. Tensor-tensor algebra for optimal representation and compression of multiway data. Proceedings of the National Academy of Sciences, 118(28):e2015851118, 2021. [27] Tamara G. Kolda and Brett W. Bader. Tensor decompositions and applications. SIAM Review, 51(3):455–500, 2009. [28] Wei Kuang, Vishwas Rao, Alexis Montoison, François Pacaud, and Mihai Anitescu. Recovering sparse DFT from missing signals via interior point method on GPU. arXiv preprint arXiv:2502.04217, 2025. [29] Canyi Lu. Tensor-Tensor Product Toolbox. Carnegie Mellon University, June 2018. https: //github.com/canyilu/tproduct. [30] Leon Mirsky. Symmetric gauge functions and unitarily invariant norms. The Quarterly Journal of Mathematics, 11(1):50–59, 1960. [31] Elizabeth Newman and Katherine Keegan. Optimal matrix-mimetic tensor algebras via variable projection. SIAM Journal on Matrix Analysis and Applications, 46(3):1764–1790, 2025. [32] Ivan V. Oseledets. Tensor-train decomposition. SIAM Journal on Scientific Computing, 33(5):2295–2317, 2011. [33] Eric T. Phipps and Tamara G. Kolda. Software for sparse tensor decomposition on emerging computing architectures. SIAM Journal on Scientific Computing, 41(3):C269–C290, 2019. GenTen: shared-memory/Kokkos parallel CP decomposition. [34] Melven Röhrig-Zöllner, Jonas Thies, and Achim Basermann. Performance of the low-rank tt-svd for large dense tensors on modern multicore cpus. SIAM Journal on Scientific Computing, 44(4):C287–C309, 2022. [35] Shaden Smith and George Karypis. A medium-grained algorithm for distributed sparse tensor factorization. In 2016 IEEE International Parallel and Distributed Processing Symposium (IPDPS), pages 902–911, 2016. SPLATT: distributed-memory parallel CP decomposition. [36] Ledyard R. Tucker. Some mathematical notes on three-mode factor analysis. Psychometrika, 31(3):279–311, 1966. [37] Sasindu Wijeratne, Rajgopal Kannan, and Viktor Prasanna. Amped: Accelerating mttkrp for billion-scale sparse tensor decomposition on multiple gpus. In Proceedings of the 54th International Conference on Parallel Processing, ICPP ’25, page 208–217, New York, NY, USA, 2025. Association for Computing Machinery.

Record · ID 194285 · SHA-256 ffdb75df79f6a36b
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.