Source-linked AI summary
Minimizing the Arithmetic and Communication Complexity of Jacobi's Method for Eigenvalues and Singular Values: Part Two -- Parallel Algorithms
James Demmel, Hengrui Luo, Ryan Schneider, Yifu Wang
TL;DR
Parallel Jacobi seeks arithmetic and communication costs comparable to distributed-memory matrix multiplication for symmetric eigenvalue and SVD problems. The paper analyzes blocked Jacobi on 2D and 2.5D processor grids, including replicated-memory communication reductions and lower bounds. It obtains perfect arithmetic scaling in 2D, near-optimal communication there, separate bandwidth- and latency-optimal 2.5D configurations, and a lower bound ruling out simultaneous best-known 2.5D matrix-multiplication communication costs.
Problem
The paper addresses how to minimize arithmetic, bandwidth, and latency costs when parallelizing Jacobi’s method for symmetric eigenvalues and SVD.
Method
The authors analyze blocked Jacobi in distributed memory on variable 2D and 2.5D processor grids, using matrix-multiplication communication strategies and proving a Jacobi-specific lower bound.
Results
2D parallel Jacobi achieves O(n^3/P) arithmetic cost, optimal bandwidth, and latency within a log P factor of the matrix-multiplication bound; 2.5D variants offer separate bandwidth- or latency-optimal configurations.
Takeaways & Limitations
Parallel Jacobi is nearly optimal in 2D, while extra memory in 2.5D can reduce one communication metric but cannot make both communication metrics simultaneously match the best-known matrix-multiplication bounds.
Takeaways & Limitations
The presented algorithms omit pivoting and therefore have no convergence guarantee; the complexity bounds are stated for a single sweep.
Abstract
from arXiv · showhide
This paper presents several parallel versions of Jacobi's method for the symmetric eigenvalue problem and the SVD. A continuation of [Demmel, Luo, Schneider, & Wang 2025], we develop parallel Jacobi algorithms whose arithmetic cost is optimal and whose bandwidth or latency can match the corresponding lower bounds of parallel matrix multiplication. Our focus is a standard distributed-memory setting with variable processor layouts, including both 2D and 2.5D processor grids. In the 2D case, we demonstrate that a standard implementation of parallel Jacobi achieves a perfect speedup in arithmetic cost -- i.e., complexity $O(n^3/P)$ when done with $P$ processors -- while hitting the 2D matrix-multiplication lower bound for bandwidth and (nearly) the lower bound for latency. By employing a 2.5D processor grid and leveraging 2.5D matrix multiplication, equivalently by increasing the memory per processor, we demonstrate that parallel Jacobi can achieve even lower bandwidth/latency, though we also prove that these costs cannot simultaneously match the best-known bounds for parallel matrix multiplication in any Jacobi algorithm. Finally, we extend our results to one-sided Jacobi SVD.
1 Introduction
The paper develops parallel Jacobi algorithms for symmetric eigenvalues and SVD, targeting arithmetic and communication costs comparable to distributed-memory matrix multiplication. It analyzes 2D and 2.5D processor layouts, proves a bandwidth–latency tradeoff, and extends the results to one-sided Jacobi SVD.
- Scope: Jacobi’s method diagonalizes symmetric matrices through orthogonal similarity transformations and can compute SVD information by diagonalizing G^T G.Its diagonal output approximates eigenvalues or squared singular values, while the accumulated orthogonal factors approximate eigenvectors or right singular vectors.
- Scope: The paper measures optimal parallel Jacobi complexity against distributed-memory parallel matrix multiplication on the same processor grid.The comparison covers arithmetic, bandwidth, and latency costs in 2D and 2.5D layouts.
- Contributions: 2D parallel Jacobi achieves arithmetic cost O(n^3/P), bandwidth O(n^2/√P), and latency O(√P log P) per sweep.The arithmetic and bandwidth costs are optimal, while latency is within a factor of log P of the 2D matrix-multiplication bound.
- Contributions: 2.5D parallel Jacobi uses replicated data and extra memory to provide separate latency-optimal and bandwidth-optimal configurations.Each processor requires at least Ω(cn^2/P) memory when the replication factor is c.
- Contributions: Any parallel Jacobi algorithm satisfying the paper’s conditions has bandwidth–latency product W(n, P)·S(n, P)=Ω(n^2), preventing simultaneous attainment of the best-known 2.5D matrix-multiplication bounds when c > 1.The paper reports that all considered Jacobi variants attain this product lower bound up to logarithmic factors.
- Contributions: The analysis is extended to parallel one-sided Jacobi SVD, while some competing eigensolver communication bounds remain open problems.The paper identifies unresolved bounds for Yau–Lu and QDWH-eig algorithms.
2 Preliminaries
The preliminaries define blocked Jacobi, its sweep-based cost model, convergence assumptions, and the distributed-memory setting used for parallelization. Blocking replaces scalar annihilations with transformations of paired off-diagonal blocks, whose asymptotic per-sweep costs are separated from the number of sweeps.
- Blocked Jacobi: Blocked Jacobi annihilates two b × b off-diagonal blocks by diagonalizing a 2b × 2b submatrix and applying the resulting orthogonal transformation to A and optionally Q.The input is partitioned into blocks AIJ with block size b dividing n.
- Cost model: A sweep passes over all off-diagonal blocks, and total cost is modeled through the number of sweeps multiplied by per-sweep messages, words moved, and floating-point operations.The sweep quantities are Ssweep, Wsweep, and Fsweep, while N is the number of sweeps needed to meet a tolerance when convergence occurs.
- Convergence: Row/column cyclic ordering and a pivoting step using QRCP or LUPP can guarantee convergence without changing the asymptotic per-sweep costs.The pivoting update applies a permutation to the local eigenvector factor.
- Convergence: The paper omits pivoting to simplify its algorithms, so its stated bounds concern one sweep and do not include a convergence guarantee.The authors report that pivoting can be added without altering single-sweep complexities and that naive implementations typically converge in O(1) sweeps.
- Distributed setting: The distributed-memory analysis uses collective operations such as scatter, broadcast, gather, and reduce on processor grids.Their associated communication costs are expressed in the same α-β-γ model.
3 2D Parallel Jacobi
The 2D algorithm distributes matrix blocks over a two-dimensional processor grid and groups independent block rotations so they execute concurrently. Its per-sweep arithmetic and bandwidth costs are optimal relative to matrix multiplication, while latency is nearly optimal.
- Processor layout: 2D parallel Jacobi assigns blocks of A and Q to a two-dimensional processor grid with b = n/√P.Each processor Proc(I, J) owns the corresponding block AIJ.
- Parallel ordering: Independent off-diagonal blocks in distinct block rows and columns can be transformed simultaneously, enabling concurrent block Jacobi updates.The algorithm groups pairwise disjoint index pairs and processes the groups cyclically by anti-diagonals.
- Algorithm: Algorithm 2 executes one sweep by computing and applying the rotations for each disjoint group in parallel, with communication updating A and the optional eigenvector matrix Q.A full computation repeatedly invokes the single-sweep routine while monitoring convergence.
- Complexity: The arithmetic speedup is perfectly linear, bandwidth reaches the matrix-multiplication lower bound under minimum memory, and latency is optimal up to log P.The algorithm is therefore characterized as nearly optimal in the 2D setting.
4 2.5D Parallel Jacobi
The 2.5D extension replicates data across processor layers and distributes Jacobi work to reduce communication. It improves either bandwidth or latency over 2D Jacobi, but not both simultaneously.
- 2.5D processor grids: 2.5D matrix multiplication lowers communication by replicating data across processor layers, at the cost of additional memory per processor.The replication factor c requires memory Ω(cn^2/P) per processor.
- Latency and bandwidth optimality: 2.5D Jacobi can improve both bandwidth and latency relative to 2D Jacobi, but its bandwidth- and latency-improvement regions are disjoint.Thus, parameter choices that optimize one communication metric increase or leave unchanged the other.
- 2.5D processor grids: 2.5D Jacobi runs the 2D algorithm on the top layer while distributing subproblem diagonalization and block updates across other layers.The input matrices A and Q are initially distributed on the top layer and replicated as needed.
- Algorithmic structure: 2.5D parallel Jacobi uses coarse B × B tiles, whose associated matrix multiplications run over μ × μ × c processor subsets, requiring c ≤ μ.Fine tiles are owned by individual processors, while coarse tiles are owned collectively by μ^2 processors.
- Communication operations: The subproblem diagonalization requires scattering each 2B × 2B subproblem across selected processors and gathering the result back to the top layer.The participating processors and block sizes depend on the utilization factor η.
- Latency and bandwidth optimality: The bandwidth-optimal configuration attains the 2.5D matrix-multiplication bandwidth bound, while latency-optimal and bandwidth-optimal configurations differ.The bandwidth-optimal choice incurs increased latency.
5 Lower Bounds for Parallel Jacobi
This section proves lower bounds for one sweep of parallel Jacobi under mild assumptions about direct diagonalization and subsequent communication. The bounds formalize unavoidable tradeoffs among arithmetic, bandwidth, and latency.
- Assumptions: The lower bounds assume direct subproblem diagonalizations occur in one processor’s fast memory and each diagonalization communicates with another processor before updates.Both assumptions are stated as mild and are satisfied by the paper’s 2D and 2.5D Jacobi algorithms.
- Implications: The resulting latency–bandwidth and latency–arithmetic tradeoffs add parallel Jacobi to numerical linear algebra computations with established tradeoff bounds.The paper reports these relationships as a broader consequence of Theorem 4.
- Lower bounds: Any qualifying parallel Jacobi algorithm satisfies W(n, P) · S(n, P) = Ω(n^2) and F(n, P) · S(n, P)^2 = Ω(n^3).Here F, W, and S denote arithmetic, bandwidth, and latency costs for one sweep.
- Proof strategy: The proof constructs a serial critical path by partitioning first-row entries into sets annihilated after successive single-processor diagonalizations.Only one collection of entries in the first row can be zeroed out at a time.
6 2.5D Parallel Jacobi SVD
The paper extends block Jacobi to a one-sided SVD algorithm that works implicitly with G^T G and parallelizes its column operations on a 2.5D grid. The construction reuses the 2.5D Jacobi machinery while avoiding explicit formation of G^T G when possible.
- Algorithm: One-sided Jacobi SVD applies block Jacobi to G^T G implicitly, forming submatrices from pairs of columns of G before updating G and V.The updates force the corresponding columns of G to become mutually orthogonal.
- Parallelization: Independent block-column pairs can be processed in parallel when their index sets are disjoint, using the usual anti-diagonal sweep ordering.This is the SVD analogue of the block-independence criterion used for parallel Jacobi.
- Processor layout: The 2.5D SVD distributes G and V across a replicated processor grid and uses coarse B × B tiles aggregated from fine b × b tiles.The coarse-tile size is B = μb, where μ is a batching parameter.
- Numerical consideration: The algorithm generally avoids explicitly forming G^T G because its condition number is the square of G’s condition number.The paper states κ_2(G^T G) = κ_2(G)^2.
- Algorithm: The algorithm targets an approximate reduced decomposition G = UΣV^T, with diagonal Σ containing approximate singular values.The columns of G converge toward UΣ while V stores approximate right singular vectors.
- Subproblem construction: Each submatrix of G^T G is computed from two coarse columns through 2.5D matrix multiplications and reductions before Jacobi diagonalization.The resulting subproblem is placed on the top-layer processor slab, after which the 2.5D Jacobi procedure continues.
I GJ and GT
The 2.5D one-sided Jacobi SVD forms small Gram subproblems with 2.5D matrix multiplication, diagonalizes them using 2D parallel Jacobi, and updates G and V. Its complexity is parameterized by batching μ and utilization η, with a latency–bandwidth tradeoff and a conjectured communication limitation.
- Complexity: Theorem 5 gives a general one-sweep complexity bound covering any batching parameter μ and utilization factor η.The analysis assumes O(1) sweeps of 2D parallel Jacobi per subproblem and accounts for the total number of parallel steps.
- Algorithmic structure: One sweep of 2.5D Jacobi SVD forms 2B × 2B Gram subproblems from paired column blocks and applies 2D parallel Jacobi to each.The subproblems use  = G(:, [I, J])^T G(:, [I, J]) and are distributed to selected processors before diagonalization.
- Algorithmic structure: The algorithm uses 2.5D matrix multiplication over μ × μ × c subgrids to form the block Gram products and then reduces partial products.Each processor participates in at most four such multiplications for a paired block-column index set.
- Complexity: Choosing fewer processors for each subproblem yields a latency-optimal implementation, while minimizing latency increases bandwidth and arithmetic costs.The paper also identifies η choices that preserve optimal arithmetic cost and discusses the analogous bandwidth–latency tradeoff for one-sided Jacobi SVD.
- Communication tradeoff: The first latency and bandwidth terms multiply to n^2 up to logarithmic factors for every μ and η, motivating conjectured analogous lower bounds for one-sided Jacobi SVD.This tradeoff persists across the parameter choices considered in Theorem 5.
7 Conclusion
The paper develops 2D and 2.5D parallel Jacobi algorithms for symmetric eigenproblems and extends the analysis to one-sided Jacobi SVD. It establishes optimal arithmetic costs and matrix-multiplication-level communication results, while identifying implementation and analysis directions that remain open.
- Conclusion: 2D and 2.5D processor grids yield parallel Jacobi algorithms with optimal arithmetic cost and bandwidth or latency matching distributed-memory matrix-multiplication bounds.The paper reports a theoretical bandwidth–latency tradeoff across these implementations.
- Conclusion: The work extends the corresponding complexity analysis to one-sided Jacobi SVD.The SVD results are presented as analogous to the symmetric eigenvalue results.
- Conclusion: The paper adds Jacobi’s method to dense parallel linear-algebra routines whose communication benefits from additional memory have been analyzed rigorously.The authors connect this result with their earlier work to reassert Jacobi’s relevance for modern scientific computing.
- Open questions: The theoretical efficiency gains of 2.5D parallel Jacobi have not yet been demonstrated in practice.The authors identify implementation, precision, and automated parameter-tuning questions as future directions.
A Detailed Pseudocode for 2D Parallel Jacobi
The appendix presents processor-level pseudocode for 2D parallel Jacobi and organizes each processor’s role through seven communication and update cases. An auxiliary bookkeeping array records the applicable cases and communication partners.
- Algorithm 8: Algorithm 8 gives a processor-level implementation of 2D parallel Jacobi using the notation introduced in Section 3.The pseudocode specifies which processors in the 2D grid each processor must communicate with.
- Processor roles: Seven cases classify whether a processor diagonalizes a block, broadcasts rotation data, or applies an update in a related row or column.The cases include processors holding blocks to be zeroed and processors involved in left-column, right-column, upper-row, or lower-row updates.
- Processor roles: Each processor stores case numbers and communication indices in an info bookkeeping array.The appendix illustrates these labels on an 8 × 8 processor grid, where participating superdiagonal processors are highlighted.
B Optimal Parameter Selection in 2.5D Parallel Jacobi
This section characterizes parameter selection and communication optimization for 2.5D parallel Jacobi. It distinguishes additive runtime from a minimax communication criterion and notes that exact machine-level optimization requires calibrated constants.
- Parameter feasibility: The feasible parameter set from Theorem 3 constrains the choices used in 2.5D parallel Jacobi.The assumption c ≥ 2 ensures that the allowable interval for η remains nonempty over the stated μ range.
- Optimization criteria: The minimax communication criterion for fixed hardware parameters α and β is distinct from the additive runtime model αS + βW + γF.The minimax formulation treats communication as a bottleneck rather than summing communication and arithmetic terms.
- Optimization criteria: Theorem 3’s big-O bounds do not determine the constant factors needed for an exact machine-level optimum.With calibrated constants C_S and C_W, the leading-term model replaces α and β by C_Sα and C_Wβ.
- Algorithm interface: The processor-level algorithm applies one sweep of block Jacobi to a distributed symmetric matrix and optionally updates an approximate eigenvector matrix Q.Its input includes A, optional Q, and a list of processor-distributed index groups.
D11 D12 D21 D22
The paper optimizes 2.5D parallel Jacobi’s latency and bandwidth through batching and utilization choices, while characterizing how the optimum changes with communication weights. The analysis also identifies practical boundaries where theoretical choices should be benchmarked or constrained.
- Parameter selection: As bandwidth becomes more important, the optimum reaches µ = c and then follows that boundary while utilization η increases.For small bandwidth weight, the optimizer may remain inside the feasible region; increasing bandwidth weight activates the lower batching constraint.
- Practical limitations: The continuous optimum is a guide rather than a prescription because hidden constants, calibration, rounding, and divisibility constraints can shift implemented choices.The recommended practice is to benchmark a small set of nearby admissible candidates on the target machine.
- Parameter selection: The optimizer is determined by the unique crossing of increasing latency and decreasing bandwidth terms, unless the crossing is clipped to a feasible endpoint.The optimization proof establishes monotonicity of the reduced terms and shows that only the ratio β/α matters.