# Multi-GPU Implementation of the Levy and Lindenbaum's Method for Streaming SVD with Application to Large Datasets of KHRTI Simulation

Yassin Ajanif, Lucas Lestandi, Luisa Silva, Aditi Sengupta, Bhavna Joshi,
Hugues Digonnet

► To cite this version:

Yassin Ajanif, Lucas Lestandi, Luisa Silva, Aditi Sengupta, Bhavna Joshi, et al. Multi-GPU Implementation of the Levy and Lindenbaum's Method for Streaming SVD with Application to Large Datasets of KHRTI Simulation. 2025. hal-05315804

HAL Id: hal-05315804

https://hal.science/hal-05315804v1

Preprint submitted on 15 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Multi-GPU Implementation of the Levy and Lindenbaum's Method for
Streaming SVD with Application to Large Datasets of KHRTI Simulation

Yassin Ajanifª,*, Lucas Lestandiª, Luisa Silvaª, Aditi Senguptaᵇ, Bhavna Joshiᵇ, Hugues Digonnetª

ªÉcole Centrale Nantes, Nantes Université, CNRS, GeM, UMR 6183, F-44000 Nantes, France

ᵇDepartment of Mechanical Engineering, Indian Institute of Technology (Indian School of Mines), Dhanbad, Jharkhand
826004, India

**Abstract**

Singular Value Decomposition (SVD) is a fundamental technique in data-driven approaches, particularly in fluid mechanics. It enables the decomposition of complex datasets, the extraction of dominant flow structures, and the construction of reduced-order models, thereby significantly reducing the complexity of physical systems. However, when applied to extremely large datasets involving millions of degrees of freedom, traditional SVD methods, even in parallel settings, become computationally prohibitive or infeasible. To overcome this limitation, the innovative approach proposed by Levy and Lindenbaum for computing low-rank approximations can be employed. This technique allows for streaming data processing, significantly reducing computational complexity in both space and time, and enabling the treatment of datasets that far exceed available memory. We propose a new parallel implementation of this method, leveraging multiple GPUs. As a case study, we apply our implementation to data obtained by Direct numerical simulation (DNS) of the Kelvin-Helmholtz-Rayleigh-Taylor instability (KHRTI), which comprises 480 million spatial degrees of freedom and 875 snapshots, demonstrating the scalability, robustness, and practical efficiency of our implementation in handling very large datasets.

Keywords:

Singular value decomposition

Proper orthogonal decomposition

High performance computing

Fluid dynamics

Large-scale simulations

This is an open access article under the **CC BY** license (http://creativecommons.org/licenses/by/
4.0/).

<footnote>*Email address: lucas.lestandi@ec-nantes.fr</footnote>
# Program Summary

**Program title:** HLL-SVD

**Developer's repository link:** https://github.com/YassinAJANIF/HLL-SVD

**Licensing provisions:** MIT license

**Programming language:** Python

**Supplementary material:** Tutorial, example datasets.

**Nature of problem:** HLL-SVD is an open-source software designed for the computation of the SVD in a streaming fashion, which is particularly advantageous when the dataset exceeds the available system memory. The software introduces a new parallel implementation of the Levy and Lindenbaum method that leverages modern GPUs, making it highly efficient for processing very large datasets, especially those arising from fluid-mechanics simulations.

**Solution method:** The HLL-SVD methodological framework is structured into two main modules. The first module, *serial_svd*, provides a serial implementation that enable the computation of the SVD in a streaming fashion. The second module, *parallel_svd*, offers parallel implementations, including our GPU-based HLL-SVD and a pure MPI version of the Levy and Lindenbaum method. It also provides additional MPI-based implementations for computing the SVD directly, without streaming. This modular design provides users with the flexibility to select the most suitable implementation according to the problem size and the available hardware resources.

**Additional comments including Restrictions and Unusual features:** HLL-SVD requires CUDA-capable GPUs. Instructions for installing dependencies and running the software are provided on the GitHub page.

# 1. Introduction

The Kelvin-Helmholtz Rayleigh-Taylor Instability (KHRTI)[1] illustrates a complex flow configuration where the interaction between shear-driven and buoyancy-driven mechanisms gives rise to rich, nonlinear, and multiscale dynamics. High-fidelity simulations of KHRTI produce massive datasets, which pose a significant challenge for both data storage and post-processing, particularly when attempting to extract the dominant mechanisms governing the system's evolution. To address this, data-driven techniques are increasingly employed to reduce the complexity of such physical systems. These approaches enable the identification of coherent structures that underlie the flow dynamics, facilitating both physical interpretation and efficient modeling. They are particularly valuable in the development of Reduced-Order Models (ROMs), which aim to reproduce the essential features of the full system at a significantly reduced computational cost.
Among the most widely adopted data-driven techniques is Proper Orthogonal Decomposition (POD) [2], which is a linear method allowing the extraction of coherent structures ranked by their $L_2$ norm, provides a compact yet accurate representation of complex flow fields. The POD also provided remarkable effectiveness in highly non-linear Dynamic system, owing to the rapid decay of the Kolmogorov N-width [3].

A direct method for computing the POD basis is the Singular Value Decomposition (SVD) [4], which operates on a snapshot matrix constructed from high-fidelity simulation data. For the latter, the Reynolds number [5] can reach values as high as $Re = 10^5$, requiring on the order of $O(Re^{9/4})$ grid points. In such cases, computing the full SVD of the snapshot matrix using classical algorithms becomes computationally infeasible. Therefore, the Truncated SVD (T-SVD) is often preferred, as it enables the extraction of only the most relevant features of the system — those that govern the dominant patterns and the dynamic behavior of the model. Given that numerical simulations are performed on powerful high-performance computing (HPC) systems and generate massive amounts of data, often hundreds of gigabytes, it becomes essential to design parallel algorithms for SVD (or T-SVD) that can fully exploit modern HPC architectures.

Sayadi et al. [6] proposed a parallel SVD implementation using the direct tall-and-skinny QR decomposition (TSQR) [7] for the data-driven decomposition of transitional boundary layers. Wang et al. [8] provided a distributed approach based on the snapshot method to approximate the POD basis. The same method was employed in the work of Farcas et al. [9] for model reduction applied to a large-scale real-world combustion scenario. In the work of Vilas [10], a parallel version of the snapshot method was presented. The author also introduced two parallel approaches to compute the Dynamic Mode Decomposition (DMD) [11, 12, 13]: the first is based on a direct TSQR algorithm, and the second on a modified TSQR algorithm [14] designed for efficient memory utilization. In the work of Iwen and Ong [15], a parallel hierarchical approach was introduced to efficiently compute the singular values and left singular vectors.

On the software implementation domain, several libraries offer parallel SVD capabilities. For example, PyLOM [16] provides an implementation of POD based on a binary tree TSQR algorithm, as in Demmel et al. [17]. The same implementation of POD was used in [18] to analyze the vortex structures in the wake of a vibrating cylinder. PyLOM includes also modules for DMD and Spectral Proper Orthogonal Decomposition (SPOD) [19], [20], and it is based on compiled functions coded with C and Cython [21] that can be called in Python scripts. Another library is libROM [22] which is a collection of C++ classes that compute reduced order models including parallel implementations like POD and DMD, another python library solution is PyDSVD-GPU [23], which computes the first $k$ monotonically decreasing singular values and the corresponding vectors using the power method [24] on multiple GPUs. Many parallel SVD algorithms require the entire dataset to be loaded into memory — whether across multiple CPU cores or GPUs. However,
this approach becomes impractical when working with extremely large datasets (e.g., hundreds of terabytes),
as is often the case in fluid dynamics applications. In such scenarios, a streaming-based approach offers
an attractive alternative by eliminating the need to load the full dataset into memory. In their seminal
paper, Levy and Lindenbaum [25](1998) proposed an innovative approach based on streaming or batched
data processing. Instead of requiring full data loading, their method processes data incrementally in smaller
batches. This approach was used by Schwerdtner et al. [26] for the construction of quadratic manifolds
from streaming data. However, when handling very large data, the serial approach can be inefficient, this
justifies the importance of parallel implementation which allows the distribution of the data on multiple cores.
Maulik and Mengaldo have developed PyParSVD [27], a python open source library based on Message Passing
Interface (MPI). This library provides a parallel implementation of the Levy and Lindenbaum approach for
computing low-rank approximations. It employs both the APMOS and direct TSQR algorithms to compute
the SVD in parallel. The library also includes a randomized SVD method [28], used to compute the SVD in
serial.

In this work, we present the Hybrid Levy-Lindenbaum SVD (HLL-SVD), a new implementation of the
Levy-Lindenbaum method designed to leverage multiple GPUs. This implementation is highly efficient in
handling large-scale datasets that exceed the system memory, such as those arising from fluid mechanics
simulations.

The proposed library, built upon the PyParSVD framework, is open source and provides various alterna-
tive parallel SVD implementations (non-streaming), such as the snapshot method, TSQR-based SVD, and
randomized SVD via TSQR. Additionally, several serial (single-core) variants are included, although they are
not within the scope of the present study.

The article is organized as follows: Section 2 describes the details of our implementation. Section 3
presents the performance results and compares the proposed approach with the previous implementation.
Finally, in Section 4, we apply our method to large-scale datasets from KHRTI simulations in order to
extract the dominant flow features.

## 2. A new implementation of the Levy and Lindenbaum approach on hybrid architectures

### 2.1. *The Levy and Lindenbaum approach*

Our primary objective in this study is to extract and analyze the dominant features of fluid mechanics
data using a streaming SVD approach, based on the method proposed by Levy and Lindenbaum. In the
classical SVD framework, the decomposition operates on a matrix $A \in \mathbb{R}^{m \times n}$, where $m$ is the number of
degrees of freedom, and $n$ is the number of snapshots. In the following, we assume that $m \gg n$, as is common
with large datasets resulting from fluid mechanics computations. The SVD factorizes the matrix A into the product of three matrices:

$$A = U \Sigma V^T \quad (1)$$

where $U \in \mathbb{R}^{m \times n}$ is the matrix of left singular vectors, $\Sigma \in \mathbb{R}^{n \times n}$ is a diagonal matrix containing the singular values in decreasing order, and $V \in \mathbb{R}^{n \times n}$ is the matrix of right singular vectors.

Algorithm 1 illustrates the functioning of the serial version of the Levy and Lindenbaum approach. The original snapshot matrix must be divided into sub-matrices (batches) with the same number of DOF and fewer columns $p \ll n$. Then, A can be written as $A = [A_1, A_2, ..., A_b]$, where $A_i \in \mathbb{R}^{m \times p}$ and $b$ is the number of batches. The value of $k$, representing the number of modes to retain, must be predefined.

The algorithm consists of two main steps. The first step is the initialization, in which the SVD of the first batch $A_1$ is computed as $A_1 = U_1 \Sigma_1 V_1^T$. This is performed by first applying a QR decomposition on $A_1 = Q_1 R_1$, followed by the SVD of the resulting matrix $R_1$. The first $k$ left singular vectors of $U_1$ and and their corresponding singular values in $\Sigma_1$ are retained, while the matrix $V_1^T$ is discarded.

The second step of the algorithm is the concatenation phase, which is then repeated iteratively for each subsequent batch. At each iteration $i$, the following matrix is constructed:

$$[ff \cdot U_{i-1} \Sigma_{i-1} | A_i]$$

where $U_{i-1} \Sigma_{i-1}$ denotes the product of the retained singular vectors and singular values from the previous iteration. A QR decomposition is applied to this concatenated matrix, followed by an SVD of the resulting upper triangular factor. Again, only the first $k$ left singular vectors $U_i$ and corresponding singular values $\Sigma_i$ are preserved for the next step.

A scalar forgetting factor $ff \in (0, 1]$ is introduced to emphasize recent data by scaling down the contribution of previously retained information. In the remainder of this paper, we fix $ff = 1$, giving equal weight to past and current data.

This approach enables the processing of large dataset without the need to retain previously seen data, thereby reducing the space complexity to $O(mk)$ and the time complexity to $O(mnk)$. In contrast, the classical approach requires $O(mn)$ memory and $O(mn^2)$ operations, making this streaming strategy significantly more efficient for large-scale datasets that exceed the memory capacity of the system.
**Algorithm 1** Streaming singular value decomposition [25]

1: **Output:** Truncated left singular vector $\mathbf{U}_i$ after $i$ batch iterations.

2: **Inputs:** list of batches $A_1, A_2, ..., A_p$.

3: **Parameters:** Forget factor $ff$;

4: **Initialization:**

5: Initial data matrix $\mathbf{A}_1 \in \mathbb{R}^{m \times b}$, where $b$ is the number of snapshots per batch;

6: Perform QR-decomposition: $\mathbf{A}_1 = \mathbf{Q}_1 \mathbf{R}_1$;

7: Perform SVD of $\mathbf{R}_1 = \mathbf{U}'\mathbf{\Sigma}_1 \mathbf{V}_1^T$ and obtain $\mathbf{U}_1 = \mathbf{Q}\mathbf{U}'$;

8: **while** New data $\mathbf{A}_i$ available **do**

9: Compute QR decomposition after column-wise concatenation of new data:

$$[ff \cdot \mathbf{U}_{i-1}\mathbf{\Sigma}_{i-1} | \mathbf{A}_i] = \mathbf{U}_{i-1}'\mathbf{\Sigma}_{i-1}'$$

10: Compute SVD of $\mathbf{\Sigma}'_{i-1} = \tilde{\mathbf{U}}_{i-1}\tilde{\mathbf{\Sigma}}_{i-1}\tilde{\mathbf{V}}^T$;

11: Preserve the first $K$ columns of $\tilde{\mathbf{U}}_{i-1}$ and denote $\tilde{\mathbf{U}}_{i-1}$;

12: Obtain the updated left singular vectors: $\mathbf{U}_i = \mathbf{U}_{i-1}'\tilde{\mathbf{U}}_{i-1}$;

13: Truncate to retain $K$ values of $\tilde{\mathbf{\Sigma}}_{i-1}$ to obtain $\mathbf{\Sigma}_i$;

14: **end while**

As mentioned in the introduction, Algorithm 1 was originally implemented by R. Maulik and G. Mengaldo [27] within an MPI-based framework that combines APMOS algorithm for the initialization phase and a direct TSQR algorithm for the concatenation step.

A major limitation of the direct TSQR approach arises in the second stage, where the local $R$ factors obtained from distributed QR decompositions are gathered onto a single processor. The size of the resulting matrix increases linearly with the number of processors, which can lead to memory bottlenecks and significant communication overhead, particularly when the local $R$ matrices are large [7]. Despite this issue, Sayadi and Schmid [6] demonstrated good scalability of the algorithm in their work, where they applied the direct TSQR method to compute POD and DMD on large-scale datasets. In Algorithm 1, the direct TSQR decomposition is used to compute the SVD in parallel during both the initialization and concatenation steps. Typically, the input matrices (batches) have a small number of columns compared to the original dataset, thus preserving the parallel efficiency of the direct TSQR algorithm.
In recent years, the emergence of new graphics cards such as the NVIDIA H100 GPUs, characterized by their high computational power and large memory bandwidth enabling fast data transfer between host and device, has made it possible to port the direct TSQR-based algorithm to these architectures. In this study, we propose to replace traditionally CPU-based computations with GPU acceleration in order to perform algebraic operations more efficiently.

### 2.1.1. HLL-SVD: Hybrid Levy and Lindenbaum SVD

To take advantage of modern software architectures and hardware accelerators, the algorithm has been implemented in Python. However, it should be noted that this readability choice does not come at the expense of performance or maintainability. Indeed, by using well maintained cupy [29] package, we can access NVIDIA's cuBLAS library and use essential linear algebra operations commonly found in LAPACK, such as EVD, or SVD directly on GPU hardware.

To enable multi-GPU computing, we use MPI bindings provided by the mpi4py library [30] as protocol of communication. Figure 1 illustrates the computation of the SVD across multiple GPUs using the direct TSQR method (Lines 7 and 8 of Algorithm 1). The computational steps are as follows:

* **Step 1:** Each GPU computes the QR decomposition of its corresponding block from the batch $A_i$, and the resulting $R_i$ matrices are gathered to processor 0.

* **Step 2:** GPU 1 performs the QR decomposition of the gathered matrix, followed by the computation of its SVD, and proceeds to Bcast the matrices $U$ and $\Sigma$ to other GPUs.

* **Step 3:** On each GPU, a matrix multiplication is performed.

* **Step 4:** In this step, the blocks $U_i$ of the left singular vector are computed.

The two matrices $U$ and $\Sigma$ are kept on the GPUs during each iteration of the algorithm, and are finally gathered onto a single MPI process on the CPU at the end.
<figure><img src="image_2.png" /><figcaption>Figure 1: SVD computation using direct TSQR on multiple GPUs.</figcaption></figure>

## 3. Performance of HLL-SVD

### 3.1. Hardware and computing environment

The first set of benchmark experiments was conducted to evaluate the performance of our implementation on a hybrid architecture. For this purpose, we used the *Jean-Zay*<fnref n="1" /> supercomputer, equipped with an *Intel Xeon Platinum 8468* processors (48 cores at 2.10 GHz), providing a total of 96 CPU cores per node, 512 GB of RAM, and 4 *NVIDIA H100* GPUs, each with 80 GB of high-bandwidth memory. *Jean-Zay* runs Red Hat version 9.4 OS and SLURM Job manager.

### 3.2. Strong Scaling Evaluation

We evaluate the performance of our implementation by taking a random matrix $A$ ($6 \times 10^7$, $3 \times 10^3$) stored in single precision (32-bit), divided into 150 batches, each consisting of 20 columns. The number of modes retained is fixed at $k = 10$. The number of GPUs ranges from 1 to 32, which is sufficient to highlight the impact of increasing parallelism on performance. Figure 2b (a) shows the execution time in seconds, and (b) shows the speed-up, computed as follows:

$$ \text{Speed-up} = \frac{T(1)}{T(N)} \qquad (2) $$

<footnote n="1">¹http://www.idris.fr/jean-zay/</footnote>
where $T(1)$ is the execution time on one GPU and $T(N)$ is the execution time on $N$ GPUs. A near-linear speedup is observed from 1 to 16 GPUs, where the efficiency is higher then 70%. At 32 GPUs, efficiency drops to around 50%. This decline is primarily attributed to the reduced workload per GPU and the increasing impact of inter-GPU communication overhead as the number of GPUs grows.

<figure><img src="image_3.png" /><figcaption>Figure 2: Strong Scaling Results</figcaption></figure>

### 3.3. Weak scaling evaluation

In this experiment, we evaluate the weak scaling performance of our implementation. We consider an initial matrix of size $(m, n) = (1 \times 10^6, 3 \times 10^3)$ stored in single precision (32-bit). We set the number of target modes to $K = 10$, and the number of columns per batch is fixed at 20. At each step of the test, the number of rows is increased proportionally to $M \times p$, where $M$ is the number of rows per GPU (kept constant), and $p$ is the number of GPUs used. This evaluation is conducted using up to 32 GPUs. An important metric in weak scaling analysis is the efficiency, defined as follows:

$$ \text{Efficiency}(p) = \frac{T_1}{T_p} \qquad (3) $$

where $T_1$ is the execution time using one GPU, and $T_p$ is the execution time using $p$ GPUs. Ideally, if the implementation exhibits perfect weak scaling, the efficiency remains close to 1 as $p$ increases. Figure 3(a) presents the execution time, while Figure 3(b) shows the corresponding efficiency. From 1 to 8 GPUs, the efficiency remains above 90%. Between 8 and 32 GPUs, it decreases gradually: approximately 89.9% for 16 GPUs and 87.3% for 32 GPUs. This slight degradation can be attributed to the growing communication overhead among the different GPUs.
<figure><img src="image_4.png" /><figcaption>Figure 3: Weak Scaling Results</figcaption></figure>

3.4. *Comparison of our implementation with the MPI version*

Table 1 presents a performance comparison between two parallel implementations of the Levy and Lindenbaum approach: PyParSVD, an existing CPU-based MPI implementation, and our custom GPU-accelerated version, HLL-SVD. The evaluation is conducted on matrices with 100 columns and a number of rows varying from $10^5$ to $10^8$, with a fixed target rank $k = 10$ and a **batch size of 10**, to enable computation on a single GPU. For small matrix sizes, PyParSVD exhibits very good performance when distributed over 48 or 96 CPU cores, often outperforming the GPU implementation due to its lower initialization overhead. However, as the matrix size increases, the GPU version becomes increasingly effective and eventually outperforms PyParSVD. Specifically, while PyParSVD takes over 2600 seconds on a single core for a matrix with $10^8$ rows, parallelization over 96 cores reduces this to 57 seconds. In comparison, our GPU version completes the same task in just 24 seconds using two GPUs. These results emphasize that our implementation is well suited for large dataset problems when sufficient GPUs are available.
<table><caption>Table 1: Performance comparison between PyParSVD and HLL-SVD on matrices (single-precision format) with 100 columns and varying numbers of rows, using *k* = 10 and a batch size of 10 columns.</caption><thead><tr><th rowspan="2">Method</th><th rowspan="2">Cores/GPUs</th><th colspan="4">Number of Rows</th></tr><tr><th>10<sup>6</sup></th><th>10<sup>7</sup></th><th>5 × 10<sup>7</sup></th><th>10<sup>8</sup></th></tr></thead><tbody><tr><td>PyParSVD</td><td>1 core</td><td>14.47</td><td>161.24</td><td>820.05</td><td>2651.38</td></tr><tr><td>PyParSVD</td><td>48 cores</td><td>0.37</td><td>10.01</td><td>57.58</td><td>113.66</td></tr><tr><td>PyParSVD</td><td>96 cores</td><td>0.24</td><td>4.74</td><td>28.44</td><td>57.10</td></tr><tr><td>HLL-SVD</td><td>1 GPU</td><td>1.38</td><td>7.83</td><td>29.06</td><td>47.45</td></tr><tr><td>HLL-SVD</td><td>2 GPUs</td><td>1.03</td><td>5.15</td><td>16.11</td><td>24.31</td></tr></tbody></table>

# 4. Application to large-scale datasets of KHRTI

## 4.1. *Description of the KHRTI Problem*

In this section, we demonstrate the efficiency of our implementation and its ability to handle extremely large datasets generated from direct numerical simulations (DNS) of the KHRTI. These datasets are issued from the study conducted by B. Joshi et al.[1], which aimed to investigate the onset and early-stage evolution of KHRTI by pressure pulse.

KHRTI emerges from the interaction between two classical hydrodynamic instabilities: the Rayleigh-Taylor instability (RTI) [31],[32], which is induced by buoyancy forces when a heavier fluid is accelerated into a lighter one, and the Kelvin-Helmholtz instability (KHI), which arises when two superposed fluid layers exhibit a velocity shear at their interface [33]. The coupling of these mechanisms significantly amplifies disturbance growth, leading to intricate interfacial dynamics and enhanced mixing—distinctive signatures of KHRTI.
<figure><img src="image_5.png" /><figcaption>Figure 4: Schematic of the computational domain for KHRTI, with labelled boundary conditions [1].</figcaption></figure>

The KHRTI simulation is carried out within a three-dimensional cuboidal domain of dimensions ($L \times H \times W$) = (0.3 m × 0.15 m × 0.1 m), as illustrated in **Fig. 3**. Initially, the domain is divided into two compartments by an insulated, planar interface. The lower compartment is filled with air at a uniform temperature of $T = 300$ K and a density of $\rho_1 = 3.6768$ kg/m$^3$, while the upper compartment contains colder air at $T = 100$ K with a density of $\rho_2 = 1.225$ kg/m$^3$. The configuration is characterized by a velocity shear across the interface: the upper fluid layer flows at a streamwise velocity of $U_1 = 3$ m/s, whereas the lower layer moves at $U_2 = 1.5$ m/s. Both fluids are assumed to behave as ideal gases. The Reynolds and Froude numbers are, respectively, $Re = 15,000$ and $Fr = 1.2365$. All the walls of the domain are adiabatically insulated, ensuring that the system remains thermodynamically isolated. Additionally, no-slip boundary conditions are applied on all solid boundaries.

At the initial time $t = 0$, a finite-thickness partition separating the two fluid streams is instantaneously displaced in the horizontal direction. This sudden removal acts as a localized perturbation, initiating the onset of mixing and triggering the development of instabilities across a wide range of spatial and temporal scales.

The simulation is based on the solution of the three-dimensional compressible Navier-Stokes equations (NSE) [34], using a computational grid with uniform spacing in the x- and z-directions, and a non-uniform grid in the y-direction to provide enhanced resolution near the top and bottom walls, as well as at the interface between the two fluid layers. Convective flux derivatives are discretized using a sixth-order non-uniform compact (NUC6) scheme [35], while boundary closures are handled using an alternate bi-diagonal ADB3 scheme [36]. Viscous fluxes are discretized with a second-order central difference scheme adapted for non-uniform grids. Time integration is performed using a four-stage, fourth-order Runge-Kutta method. The simulation employs approximately 480 million grid points and a time step of $\Delta t = 6.25 \times 10^{-7}$ s. Figure 5
shows the evolution of the pressure disturbance during the onset of the KHRTI in the (z = 10) plane, computed as in [1]. At t = 0.0003001 in Fig. 5(a), it illustrates the formation of pressure waves at the interface separating the two fluids. These waves propagate normally to the interface, in agreement with the results of A. Sengupta [37] on the onset of the RTI triggered by acoustic excitation. At t = 0.0195039 in Fig. 5(b), the propagation of the pressure wave becomes radial, emanating from the right boundary toward the left. This behavior is induced by the shearing action of the KHI mechanism, highlighting its influence on the RTI mechanism during the onset of the KHRTI. At t = 0.0387077 in Fig. 5(c), the pressure disturbance propagates throughout the computational domain. The initially radial wave fronts reflect from the upper and lower boundaries, leading to the formation of chaotic structures, particularly visible near the left side of the domain [1].

<figure><img src="image_6.png" /><figcaption>Figure 5: The evolution of pressure disturbances in z = 10 plane during early onset of the KHRTI.</figcaption></figure>
4.2. SVD Computation for the KHRTI Dataset

As previously mentioned, the primary objective of the study [1] is to investigate the early-stage dynamics
of the KHRTI by analyzing the evolution of the pressure disturbance. Using SVD to extract the POD basis
from this field offers valuable insight into the influence of shear and buoyancy-driven instabilities on the
pressure distribution by capturing their underlying spatiotemporal dynamics.

To compute the POD basis using our implementation, we utilize 875 snapshots of the pressure disturbance field, each containing approximately 4.8×10⁸ degrees of freedom. The full dataset is stored in double precision and occupies approximately 2.4 TB of storage. For the preprocessing step, the dataset is divided into batches of 25 snapshots, resulting in a total of 35 batches. A target rank of $k = 10$ is chosen to enable the extraction of the dominant structures that capture most of the system’s energy [38]. The computation proceeds in single precision (running on 8 GPUs): a first pass computes the spatial modes, followed by a second pass that extracts the temporal modes representing the time evolution of the spatial basis.

<figure><img src="image_7.png" /><figcaption>Figure 6: Singular values of the KHRTI data, computed using the snapshot method.</figcaption></figure>

Figure 6 shows the singular-value spectrum obtained with the snapshot method in which a rapid, near-
exponential drop over the first 30 modes is observed, followed by a slow decay extending to the last nonzero
value ( 703), beyond this index, the singular values are numerically zero due to machine precision. These
results are interesting for future study and motivate a POD-Galerkin reduced-order model as a promising
approach to our problem. Figure 7 presents the first six temporal modes, computed as $a_i(t) = \sigma_i v_i(t)$, where
$\sigma_i$ denotes the i-th singular value and $v_i$ the corresponding temporal mode, they describe the time evolution
of the corresponding spatial structures shown in Figure 8. To further analyze their dynamics, we also include
their frequency content in Figure 9, obtained using the Fast Fourier Transform (FFT), which reveals the dominant frequencies and the oscillatory behavior associated with each mode.

Modes 1 and 2 form a typical oscillatory pair [39], displaying a nearly symmetric sinusoidal pattern with an amplitude that increases over time. Their frequency spectra reveal a clear peak around 0.960 kHz, with no significant low-frequency components. This indicates a purely oscillatory signal, associated with pressure waves generated by KHRTI. Modes 3 and 4, in contrast, exhibit a slower temporal evolution, with very few oscillations. As can be seen in their corresponding FFT spectra, this behavior suggests the presence of a slow, non-oscillatory component; thus, these are often referred to as ‘shift modes’ [38]. Modes 5 and 6 also form a pair modes, they display transient oscillations, localized in time and followed by a progressive decay. Their frequency spectra show a broader and less distinct peak centered around 0.971 kHz, which suggests the presence of non-stationary pressure waves, possibly associated with secondary perturbations in the system.

Figure 8 presents the first six spatial modes in (x, y)-plan with $z = 10$, associated with the temporal modes shown in Figure 7. Modes 1 and 2 exhibit similar spatial structures, characterized by a radial propagation of pressure waves from the right boundary, which gradually becomes chaotic in the left part of the domain, as illustrated in Figure 5 (c) of the pressure disturbance evolution. These two modes form a pair of modes with a phase shift and capture the dominant oscillations induced by the KHRTI phenomenon. Modes 3 and 4 reveal a structure resembling the initial configuration of the two fluids: pronounced variations appear along the y-axis, while variations along the x-axis remain weaker. Finally, for modes 5 and 6. Their spatial structure, are very similar to that of modes 1 and 2, highlights the presence of secondary oscillations within the system.

Another important role of the SVD is its ability to perform dimensionality reduction. Figure 10 illustrates the reconstruction of a snapshot at time $t = 0.02050410$ using only $k = 10$ modes. The relative error of the reconstruction is computed as follows:

$$ \varepsilon_{\text{rel}} = \frac{\|A_j - \hat{A}_j\|_2}{\|A_j\|_2}, \qquad (4) $$

where $A_j$ denotes the original snapshot with index $j$ (in our case $j = 200$), and $\hat{A}_j$ represents the reconstructed snapshot. In this case, the measured relative error is $\varepsilon_{\text{rel}} = 7.5584\%$. This value indicates that the main features of the data are well preserved, as can be observed in Figure 10.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>Figure 7: Temporal evolution of the first six POD modes for the KHRTI data.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 8: Contour plots of first six spatial/eigenfunction modes from POD of KHRTI data, with z=10.</figcaption></figure>
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>Figure 9: FFT of the first six temporal modes from POD for the KHRTI data.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 10: Reconstruction of the snapshot at $t = 0.02050410$ s using $k = 10$ modes, in the $z = 10$-plane.</figcaption></figure>

Conclusion

In this work, we have presented a novel implementation of the Levy and Lindenbaum approach, specifically designed to process extremely large datasets that exceed the memory capacity of a single system. We demonstrated the scalability of our implementation through the use of a direct TSQR algorithm on modern hybrid architectures. Our implementation leverages MPI for inter-process communication and CUDA for accelerating linear algebra operations, such as, QR decomposition, SVD, and matrix multiplications on GPUs. The implementation was validated using large-scale data from a direct numerical simulation (DNS) of the KHRTI, successfully capturing the dominant flow structures associated with pressure disturbances. This provides valuable insight into the underlying dynamics of the system. As a future direction, we aim to further enhance the performance of the Levy and Lindenbaum framework by exploring alternative parallel TSQR algorithms optimized for modern GPU architectures. This would enable more efficient utilization of computational resources and improve the overall scalability of the algorithm.
Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Code availability

The codes used for this work are available at : https://github.com/YassinAJANIF/HLL-SVD.git

Data Availability

The data that support the findings of this study are available from Lucas Lestandi (lucas.lestandi@ec-nantes.fr) and Aditi Sengupta (aditi@iitism.ac.in) upon reasonable request.

Declaration of generative AI and AI-assisted technologies in the writing process

During the preparation of this work the authors used ChatGPT (OpenAI) in order to improve the read-
ability and clarity of the manuscript. After using this tool, the authors reviewed and edited the content as
needed and take full responsibility for the content of the published article.

Acknowledgements

This work was supported by the CEFIPRA (Indo-French Centre for the Promotion of Advanced Research) within the framework of an Indo-French collaborative project. It also received support from a French government grant managed by the Agence Nationale de la Recherche (ANR) under the "Investissements d'avenir" program (reference ANR-21-ESRE-0051). Access to the MesoNET resources center and the MesoNET project was granted under allocation M24080. We gratefully acknowledge the computational resources provided by the GLiCID local HPC cluster at Centrale Nantes, as well as access to the Jean Zay supercomputing facility, provided by GENCI (Grand Équipement National de Calcul Intensif) under project allocation [AD010615966].
References

[1] Bhavna Joshi, Tapan K. Sengupta, Prasannabalaji Sundaram, and Aditi Sengupta. Highly resolved peta-scale direct numerical simulations: Onset of Kelvin-Helmholtz Rayleigh-Taylor instability via pressure pulses. *Computers & Fluids*, 284:106442, November 2024.

[2] Gal Berkooz, Philip Holmes, and John Lumley. The proper orthogonal decomposition in the analysis of turbulent flows. *Annual Review of Fluid Mechanics*, 25:539-575, 2003.

[3] Andrei N. Kolmogorov. Über die beste annäherung von funktionen einer gegebenen funktionenklasse. *Annals of Mathematics*, 37(1):107-110, 1936.

[4] Anindya Chatterjee. An introduction to the proper orthogonal decomposition. *Current Science*, 78(7):808-817, 2000.

[5] William Cabot and Andrew Cook. Reynolds number effects on rayleigh-taylor instability with possible implications for type ia supernovae. *Nature Physics*, 2:562-568, 07 2006.

[6] Taraneh Sayadi and Peter J. Schmid. Parallel data-driven decomposition algorithm for large-scale datasets: with application to transitional boundary layers. *Theor. Comput. Fluid Dyn.*, 30(5):415-428, October 2016.

[7] Austin R. Benson, David F. Gleich, and James Demmel. Direct qr factorizations for tall-and-skinny matrices in mapreduce architectures. In *2013 IEEE International Conference on Big Data*, page 264-272. IEEE, October 2013.

[8] Zhu Wang, Brian McBee, and Traian Iliescu. Approximate partitioned method of snapshots for POD. *Journal of Computational and Applied Mathematics*, 307:374-384, December 2016.

[9] Ionut-Gabriel Farcaş, Rayomand P. Gundevia, Ramakanth Munipalli, and Karen E. Willcox. Distributed computing for physics-based data-driven reduced modeling at scale: Application to a rotating detonation rocket engine. *Computer Physics Communications*, 313:109619, 2025.

[10] Vilas Shinde. Distributed-parallel proper orthogonal/dynamic mode decompositions of large flow data. *Computer Physics Communications*, 313:109644, 2025.

[11] Peter Schmid. Dynamic mode decomposition of numerical and experimental data. *Journal of Fluid Mechanics*, 656.
[12] Jonathan H. Tu, Clarence W. Rowley, Dirk M. Luchtenburg, Steven L. Brunton, and J. Nathan Kutz. On dynamic mode decomposition: Theory and applications. *Journal of Computational Dynamics*, 1(2):391–421, 2014.

[13] J. Nathan Kutz, Steven Brunton, Bingni Brunton, and Joshua Proctor. *Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems*. 11 2016.

[14] Mark Hoemmen. A communication-avoiding, hybrid-parallel, rank-revealing orthogonalization method. In *Proceedings of the 2011 IEEE International Parallel & Distributed Processing Symposium*, IPDPS '11, page 966–977, USA, 2011. IEEE Computer Society.

[15] Mark Iwen and B. W. Ong. A distributed and incremental svd algorithm for agglomerative data analysis on large networks. *SIAM Journal on Matrix Analysis and Applications*, 37(4):1699–1718, January 2016.

[16] Benet Eiximeno, Arnau Miró, Beka Begiashvili, Eusebio Valero, Ivette Rodriguez, and Oriol Lehmkhul. PyLOM: A HPC open source reduced order model suite for fluid dynamics applications. *Computer Physics Communications*, 308:109459, March 2025.

[17] James Demmel, Laura Grigori, Mark Hoemmen, and Julien Langou. Communication-optimal parallel and sequential qr and lu factorizations: theory and practice, 2008.

[18] Benet Eiximeno, Arnau Miró, Juan Carlos Cajas, Oriol Lehmkuhl, and Ivette Rodriguez. On the wake dynamics of an oscillating cylinder via proper orthogonal decomposition. *Fluids*, 7(9), 2022.

[19] Aaron Towne, Oliver T. Schmidt, and Tim Colonius. Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis. *Journal of Fluid Mechanics*, 847:821–867, 2018.

[20] Gianmarco Mengaldo and Romit Maulik. Pyspod: A python package for spectral proper orthogonal decomposition (spod). *Journal of Open Source Software*, 6(60):2862, 2021.

[21] Cython Project. Cython: C-extensions for python, 2020.

[22] Youngsoo Choi, William Arrighi, Dylan Copeland, Robert Anderson, and Geoffrey Oxberry. libROM: Library for reduced order modeling. Computer Software, 2019.

[23] Ismael Boureima, Manish Bhattarai, Maksim E. Eren, Nick Solovyev, Hristo Djidjev, and Boian S. Alexandrov. Distributed out-of-memory svd on cpu/gpu architectures, 2022.

[24] Abdeslem Bentbib and A. Kanber. Block power method for svd decomposition. *Analele Stiintifice ale Universitatii Ovidius Constanta, Seria Matematica*, 23:45–58, 06 2015.
[25] Avraham Levy and Michael Lindenbaum. Sequential Karhunen-Loeve basis extraction and its application to images. In *Proceedings 1998 International Conference on Image Processing. ICIP98 (Cat. No.98CB36269)*, volume 2, pages 456–460, Chicago, IL, USA, 1998. IEEE Comput. Soc.

[26] Paul Schwerdtner, Prakash Mohan, Aleksandra Pachalieva, Julie Bessac, Daniel O'Malley, and Benjamin Peherstorfer. Online learning of quadratic manifolds from streaming data for nonlinear dimensionality reduction and nonlinear model reduction, 2024.

[27] Romit Maulik and Gianmarco Mengaldo. PyParSVD: A streaming, distributed and randomized singular-value-decomposition library, August 2021. arXiv:2108.08845 [cs].

[28] Nathan Halko, Per-Gunnar Martinsson, and Joel A. Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions, 2010.

[29] Ryosuke Okuta, Yuya Unno, Daisuke Nishino, Shohei Hido, and Crissman Loomis. CuPy: A numpy-compatible library for NVIDIA GPU calculations. In *Proceedings of the Workshop on Machine Learning Systems (LearningSys) at the 31st Annual Conference on Neural Information Processing Systems (NeurIPS)*, 2017.

[30] Lisandro Dalcin and Yao-Lung Leo Fang. mpi4py: Status update after 12 years of development. *Computing in Science and Engineering*, 23(4):47–54, 2021.

[31] Lord Rayleigh. Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density. *Proceedings of the London Mathematical Society*, 14(1):170–177, 1882. Web.

[32] Geoffrey Ingram Taylor. The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. i. *Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences*, 201:192–196, 1950.

[33] Ye Zhou. *Hydrodynamic Instabilities and Turbulence: Rayleigh-Taylor, Richtmyer-Meshkov, and Kelvin-Helmholtz Mixing*. Cambridge University Press, December 2024.

[34] Klaus A. Hoffmann and Steve T. Chiang. *Computational Fluid Dynamics*, volume 1 of *Computational Fluid Dynamics*. Engineering Education System, 2000. Accessed via Google Books.

[35] Nidhi Sharma, Aditi Sengupta, Manoj Rajpoot, Roshan J. Samuel, and Tapan K. Sengupta. Hybrid sixth order spatial discretization scheme for non-uniform cartesian grids. *Computers & Fluids*, 157:208–231, 2017.
[36] Tapan K. Sengupta and Aditi Sengupta. A new alternating bi-diagonal compact scheme for non-uniform grids. *Journal of Computational Physics*, 310:1-25, 2016.

[37] Aditi Sengupta, Prasannabalaji Sundaram, Vajjala K. Suman, and Tapan K. Sengupta. Three-dimensional direct numerical simulation of Rayleigh-Taylor instability triggered by acoustic excitation. *Physics of Fluids*, 34(5):054108, May 2022.

[38] Bernd Noack, Konstantin Afanasiev, Marek Morzynski, Gilead Tadmor, and Frank Thiele. A hierarchy of low-dimensional models for the transient and post-transient cylinder wake. *Journal of Fluid Mechanics*, 497:335-363, 2003.

[39] Tapan K. Sengupta. *Instabilities of Flows and Transition to Turbulence*. CRC Press, Boca Raton, USA, 2012.