Extracting the Potential of Emerging Hardware Accelerators for Symmetric Eigenvalue Decomposition

The tridiagonalization process is usually the pre-step to eigenvalue decomposition. The goal of tridiagonalization can be expressed as follows:

where $Q$ is an orthogonal matrix, and $T$ is a tridiagonal matrix. The conventional tridiagonalization (Dongarra et al., 1989) uses Householder reflection to eliminates the elements except the diagonal and the subdiagonal elements in one column. Unfortunately, on the modern computer architectures, the tridiagonalization is not efficient at all, because it contains too many BLAS2 operations (such as symmetric matrix vector multiplication (symv)), which are not efficient and takes over 90% of the total elapsed time. To solve this issue, 2-stage tridiagonalization is proposed and it’s has a better utilization of hardware resources.

The 2-stage tridiagonalization has two successive processes, the first is called successive band reduction (SBR) that reduces the matrix to a band form: $A=Q^{-1}\times B\times Q$, where $Q$ is an orthogonal matrix and $B$ is a band matrix with bandwidth $b$. The 2nd stage is bulge chasing, it converts the band form matrix $B$ from 1st stage to a tridiagonal matrix $T$, see Figure 1 for a intuitive illustration. The 2-stage tridiagonalization is mathematically equivalent to direct tridiagonalization, but it provides more BLAS3 operations in SBR. Although the bulge chasing has many BLAS2 operations, the time complexity is much lower than SBR, so it doesn’t take too much time to be performed. Therefore, the 2-stage is proved to be more efficient on multi-core architectures (Haidar et al., 2011; Gates et al., 2018).

After obtaining the tridiagonal matrix, iterative methods are used to get eigenvalues. The most popular iterative methods for computing eigenvalues and eigenvectors are the QR algorithm (Watkins, 1982) and divide and conquer (D&C) (Tisseur and Dongarra, 1999), both of which are included in numerical linear algebra libraries such as cuSOLVER. The primary difference between these methods is their computational complexity: if only eigenvalues are needed, the QR algorithm requires $O(n^{2})$ operations,whereas D&C requires $O(n^{3})$ because it also generates eigenvectors during the iterations but it has better parallelism than QR algorithm.

The emerging hardware accelerators are designed to handle the computational intensity of DNNs, especially Transformer-based DNNs (Vaswani et al., 2017), which typically involve a substantial amount of matrix multiplications (GEMMs). For instance, Nvidia’s latest GPUs feature Tensor Cores that excel in low or mixed-precision GEMMs, offering nearly 1 PFLOPs peak performance for FP16 precision GEMMs on H100 GPU. Additionally, high-precision (FP64) GEMMs can also be executed on Tensor Cores, delivering performance comparable to FP32 precision GEMMs.

While the computing capacity of hardware continues to increase rapidly, the memory bandwidth is also improving, but at a slower pace. This discrepancy leads to a widening memory-compute gap.

To better understand the difficulties in fully utilizing the performance capabilities of emerging hardware, we employ the roofline model (Williams et al., 2009) and present Figure 3 to illustrate the performance differences across various generations of Nvidia GPUs. As depicted, achieving peak performance on the latest GPUs, such as the H100 and A100, necessitates a higher data intensity. Having these features, however, conventional numerical linear algebra algorithms, like the EVD solver that utilize blocking algorithms, do not achieve the requisite data intensity to fully exploit the computational power of these advanced GPUs, which means compute-bound algorithms can be converted to memory-bound on new hardware. Thus, to understand the reasons behind the above conversion, we use EVD solver as an example and address the performance bottlenecks in the following section.

Extensive literature indicates that two-stage tridiagonalization outperforms direct tridiagonalization, even when forming eigenvectors (Haidar et al., 2011; Ltaief et al., 2012; Luszczek et al., 2011). Thus, we will not delve deeply into the comparative analysis of these two approaches, but rather, we will focus on the bottlenecks of the state-of-the-art 2-stage tridiagonalization algorithm.

In this section, we discuss the hardware features of the latest GPUs: the memory speed scales much slower than that of computations for evolving GPUs, which results in an enlarging gap between compute capacity and memory bandwidth. This gap necessitates higher data intensity to achieve peak performance, complicating algorithm design. Based on our analysis of the performance bottlenecks, the conventional EVD algorithms, which is efficient on previous computer architectures, fails to handle the larger memory-compute gap and rich parallelism of emerging hardware, leading to a memory-bound implementations.

As a result, the features and bottlenecks leave us some opportunities to outperform the SOTA algorithms which are memory-bound on emerging hardware accelerators. In terms of EVD solver, in cuSOVLER and MAGMA implementations, given a matrix size $50K\times 50K$ on H100 GPU, the tridiagonalization process takes over 95% (78.5 seconds out of 80.4 seconds), and 60% (43.5 seconds out of 73.5 seconds), respectively. Thus, the EVD solver’s performance will boost significantly if we can optimize the tridiagonalization quite well by utilizing the computing capacity and parallelism on new hardware. And we’ll propose a new method of tridiagonalization to increase data intensity and handle the memory-bound problem in the following section.

Each iteration in the conventional band reduction algorithm involves a tall and skinny QR (TSQR) factorization and trailing matrix updates using the syr2k routine. There is extensive literature on performing fast TSQR on GPUs (Zhang et al., 2020a; Anderson et al., 2011; Ballard et al., 2014), which we can leverage directly. With the TSQR algorithm, the critical path is the trailing matrix update. As shown in Table 1, typical bandwidth selections of 128 or 256 achieve less than 30% of peak performance on the H100 GPU. Therefore, significantly improving syr2k’s performance would yield substantial gains in band reduction.

One straightforward method is to increase the bandwidth $b$. From Table 1, only a $k$ value of 1024 or larger in the syr2k operation achieves satisfactory trailing matrix update performance. However, increasing the bandwidth indiscriminately degrades the bulge chasing performance, as seen in Figure 4 where even $b=128$ significantly increases the bulge chasing time cost.

To address this issue, we propose a new band reduction algorithm called Detached Band Reduction (DBR), and the key idea is to decouple the bandwidth $b$ from the blocksize $nb$ (Figure 5). For a matrix $A\in\mathbb{R}^{n\times n}$, with bandwidth $b$ and block size $nb$, the algorithm is illustrated in Algorithm 1. If $b$ equals $nb$, DBR degrades to SBR.

DBR’s primary advantage over SBR is the decoupling of blocksize $nb$ from bandwidth $b$, offering greater flexibility in selecting and adjusting these parameters. In SBR, the typical choice is $b=nb=64$, whereas DBR allows configurations like $b=32$ and $nb=1024$. This algorithmic optimization provides two significant performance benefits: 1) enabling much larger $k$ values in trailing matrix updates using the syr2k routine; 2) allowing smaller bandwidth $b$ to reduce the time complexity of the subsequent bulge chasing process. These advantages lead to substantial performance improvements in band reduction, even with smaller bandwidths.

It is widely accepted that the two-stage tridiagonalization is faster on modern architectures because it converts many BLAS2 operations into BLAS3 operations. However, previous research asserts that the bulge chasing process, being limited in parallelism and close to memory bandwidth limits, would not benefit significantly from an accelerator-based implementation, as it is already optimized for CPU caches (Section 4.2 in (Gates et al., 2018)). This explains why existing 2-stage tridiagonalization/bidiagonalization algorithms use CPUs for bulge chasing (Luszczek et al., 2011; Gates et al., 2018).

In this section, we refute these claims by demonstrating that the bulge chasing process, despite being memory-bound, still exhibits enough parallelism to benefit significantly from hardware accelerators. This parallelism is evident in two aspects: 1) inter-kernel parallelism (e.g., pipelining); 2) intra-kernel parallelism (e.g., Householder transformations).

The inter-kernel parallelism is enabled because different sweeps lack data dependencies, allowing for pipelining. To illustrate this pipeline, we categorize operations into three types: Householder vector generation, matrix update from both sides, and matrix update from the left side (creating a bulge), denoted by different colors in Figure 6. We observe that when the $i$-th sweep completes three cycles, the $(i+1)$-th sweep can safely start without any data dependency. This allows one thread block to handle one sweep while other thread blocks wait and execute sequentially.

The intra-kernel parallelism is primarily leveraged in Householder vector generation and Householder transformations. By launching multiple threads, the transformation can be computed in parallel to achieve satisfactory performance. See Algorithm 2 for details.

Despite being memory-bound rather than compute-bound, bulge chasing can significantly benefit from emerging hardware accelerators, such as H100, A100, and RTX 4090 GPUs, by fully exploiting the parallelism between different sweeps and within CUDA kernels. Compared to the CPU implementation, the bulge chasing using hardware accelerators is nearly 8.0x faster than the CPU-based implementation (Figure 9).

As shown in Algorithm 1, the first step is panel QR factorization. Extensive research exists on implementing TSQR on GPUs, so we will not delve into it here; see (Zhang et al., 2020a; Ballard et al., 2014) for high-performance TSQR implementations and the reconstruction of the WY representation.

The subsequent step is panel updates between $b$ and $nb$ (Line 6 in Algorithm 1). Similar to syr2k performance (Table 1), The GEMMs shapes and sizes also affect the performance.the shapes and sizes of GEMMs significantly affect performance. Unfortunately, the GEMMs in panel updates typically have a dimension equal to $b$, resulting in suboptimal performance. To address this, we can use a recursive strategy to enlarge these dimensions.

For illustration, consider $b=32$ and $nb=256$ in DBR. Non-recursive panel updates generate 7 GEMMs with $k=32$, while recursive panel updates produce 4 GEMMs with $k=32$, 2 GEMMs with $k=64$, and 1 GEMM with $k=128$. This modification does not increase the total computations but converts GEMMs into more square shapes, increasing FLOPs efficiency.

Although algorithmic design can speed up trailing matrix updates using cuBLAS Dsyr2k routines, our experimental results show that cuBLAS routines still consume excessive time (over 90% of DBR time). Therefore, optimizing syr2k is crucial.

Testing cuBLAS Dsyr2k on an H100 GPU reveals it achieves less than 60% (36 TFLOPs) of the H100 peak FP64 precision. Additionally, a severe bug reduces Dsyr2k performance to 4 TFLOPs when $n\geq 49152$.

To improve performance, we propose a recursive formulation that maximizes the size of internal GEMMs:

This results in:

The original syr2k decomposes into two sub-syr2k operations and one large GEMM. Recursive decomposition improves hardware utilization by creating larger GEMMs. However, pure recursion reduces parallelism. Therefore, we adapt the recursive idea into a parallel version, ensuring sub-syr2k operations run in parallel. Figure 7 illustrates the iterative approach for computing syr2k. In the first iteration, diagonal blocks are computed in parallel (using cuBLAS batched GEMM), followed by iterative computation of off-diagonal blocks until the largest GEMM is processed. See Algorithm 3 for details.

The speedup compared to cuBLAS Dsyr2k routine is shown in Figure 8. We compared the square SYR2K and tall and skinny SYR2K and the experiments reveal that we can outperform cuBLAS on varies sizes and shapes.

From our previous analysis and Algorithm 2, the design strategy for bulge chasing is straightforward: each kernel handles one sweep.

Within the kernel, multiple threads perform the Householder transformations from both the left and right sides, with the number of threads determined by the bandwidth $b$. Given that bulge chasing is memory-bound, it is crucial to hide the data movement between shared and global memory. We achieve this by maintaining two shared memory blocks: one for online computations and another for data movement. Once the computation block completes its tasks, its role switches to data movement, while the other block takes over computations. This simple pipelining allows for overlapping some data movements.

A critical issue between different sweeps (kernels) is synchronization due to data dependencies between successive sweeps. Specifically, the $i+1$-th sweep can only commence after the third bulge elimination in the $i$-th sweep. Nvidia’s cooperative group tool facilitates synchronization between threadblocks. However, for this to function correctly, the number of launched blocks must not exceed the device limit, or the cooperative group will return an error indicating ”too many blocks in cooperative launch cudaLaunchCooperativeKernel.”

In our implementation, to maximize parallelism, the bulge chasing process launches a large number of threadblocks, particularly for large matrices. Therefore, using a cooperative group is impractical in this scenario. To resolve the synchronization issue, we manually set up locks between adjacent sweeps to avoid conflicts. The $i$-th sweep shares a lock flag with the $i+1$-th sweep to indicate the current column being processed. After the $i+1$-th sweep eliminates a bulge, it waits until there is no conflict as indicated by the lock. Similarly, the $i+2$-th sweep monitors the lock flag shared with the $i+1$-th sweep.

Our proposed bulge chasing implementation and optimization take full advantage of hardware accelerators, as shown by the performance comparison with MAGMA in Figure 9. The proposed implementation achieves 8.0x speedup compared to the CPU-based version.

Like conventional SBR, tuning parameters such as $b$ and $nb$ is essential to achieve optimal performance on different hardware. The bandwidth $b$ is crucial for balancing bulge chasing and trailing matrix updates. In our DBR, $nb$ must also be carefully selected to ensure peak performance during trailing matrix updates. For instance, on an A100 GPU, setting $nb=512$ is sufficient for high-performance syr2k operations, and increasing $nb$ further would result in slower panel updates. Conversely, on an H100 GPU, syr2k operations require $nb\geq 2048$, so we set $nb=2048$. In general, the block size $nb=k$, where $k$ in syr2k operation yields the best performance for large matrices.

Although we have adapted panel updates using a recursive formulation, the bandwidth $b$ does not significantly impact DBR performance, as many tall and skinny GEMMs are converted into relatively square GEMMs. However, GEMMs with one small dimension still exist and are non-negligible. Additionally, bulge chasing performance is closely related to bandwidth $b$, with smaller $b$ being preferred. Thus, it is important to find a balance between DBR and bulge chasing through tuning. Fortunately, DBR is not highly sensitive to tuning, and even suboptimal settings can deliver satisfactory performance. See Table 2) for a performance reference.

In this section, we depict our implementations and optimizations of the tridiagonalization process, including DBR, BLAS3 operations, and bulge chasing. We’ve observed significant acceleration compared to MAGMA on the H100 GPU.

The final step of the EVD solver is computing the eigenvalues using iterative methods such as QR algorithm. As previously discussed, the iterative method is not a bottleneck on GPUs. For example, cuSOLVER’s well-optimized divide and conquer only takes about 3% of the time in the Dsyevd routine. We acknowledge that outperforming cuSOLVER’s implementation in this regard is unlikely. However, due to the tremendous speedup in tridiagonalization, even a sub-optimal divide and conquer implementation can still result in substantial acceleration of the full eigenvalue decomposition.

Another observation is that our proposed methods are also beneficial on emerging GPU architectures such as the RTX 4090 series, which have limited FP64 precision computing capacity (only $\frac{1}{64}$ of FP32 peak performance). Our methods can take advantage of new technologies that leverage INT8 Tensor Cores to perform FP64 GEMMs (Ootomo et al., 2024). Performance results will be shown in the experiments section.

The boost of hardware accelerators’ developments began with Nvidia’s release of the P100 GPU in 2016, which supported half-precision computations (Turner et al., 2018). Driven by the increasing demand for training larger neural networks, the Tesla V100 GPU (Martineau et al., 2018) introduced Tensor Cores (Markidis et al., 2018), enabling exceptionally fast fused GEMMs. High precision computations have similarly benefited from these innovations. On the A100 (Choquette et al., 2021) and H100 GPUs (Choquette, 2023), FP64 GEMMs can also utilize Tensor Cores, achieving performance levels equivalent to FP32 computations. Specifically, for scientific computing requiring high precision, AMD’s recently introduced MI300x (Smith et al., 2024) offers the highest FP64 computing capacity, with a peak performance almost three times that of the H100 GPU, reaching 163.4 TFLOPs.