Source-linked AI summary
Communication-optimal parallel and sequential QR and LU factorizations
James Demmel, Laura Grigori, Mark Hoemmen, Julien Langou
TL;DR
Dense QR factorization needs communication-efficient parallel and sequential algorithms, but existing LAPACK and ScaLAPACK methods communicate asymptotically more. The paper develops TSQR and CAQR, proves communication lower bounds for LU and QR, and shows its algorithms attain those bounds up to polylogarithmic factors while matching Householder QR’s stability. These results identify communication-optimal QR algorithms and relate them to lower-bound-attaining LU methods.
Problem
Existing LAPACK and ScaLAPACK QR algorithms perform asymptotically more communication than the communication lower bounds targeted for dense QR factorization.
Method
The paper develops TSQR and CAQR and extends matrix-multiplication bandwidth lower bounds to latency lower bounds that apply to dense LU and QR.
Results
The QR algorithms attain communication lower bounds up to polylogarithmic factors, with parallel TSQR using log P messages and sequential TSQR transferring 2mn words.
Takeaways & Limitations
Communication-optimal QR is achievable for tall-skinny and general rectangular matrices while retaining conventional Householder QR numerical stability.
Takeaways & Limitations
The algorithms and optimality proofs are analyzed in detail for two-level memory hierarchies or homogeneous parallel machines, not fully general architectures.
Abstract
from arXiv · showhide
We present parallel and sequential dense QR factorization algorithms that are both optimal (up to polylogarithmic factors) in the amount of communication they perform, and just as stable as Householder QR. We prove optimality by extending known lower bounds on communication bandwidth for sequential and parallel matrix multiplication to provide latency lower bounds, and show these bounds apply to the LU and QR decompositions. We not only show that our QR algorithms attain these lower bounds (up to polylogarithmic factors), but that existing LAPACK and ScaLAPACK algorithms perform asymptotically more communication. We also point out recent LU algorithms in the literature that attain at least some of these lower bounds.
1. Introduction.
The paper develops TSQR and CAQR algorithms that minimize communication costs for parallel and sequential dense QR while retaining Householder-level numerical stability. It establishes communication lower bounds and shows the new algorithms meet them, whereas LAPACK and ScaLAPACK generally communicate more.
- Algorithms: TSQR targets tall-skinny matrices, while CAQR targets general rectangular matrices using 1-D and 2-D layouts, respectively.TSQR also serves as a building block for CAQR.
- Algorithms: The algorithms are as numerically stable as conventional Householder QR and have the same highest-order floating-point operation terms as LAPACK and ScaLAPACK.They may perform slightly more floating-point operations overall.
- Parallel communication: log P messages make parallel TSQR optimal and a factor 2n fewer than ScaLAPACK’s parallel QR factorization.This comparison assumes m/P ≥ n.
- Parallel communication: Θ(nP/m) messages make parallel CAQR optimal and a factor Θ(mn/P) fewer than ScaLAPACK on general rectangular matrices.The reported message count ignores polylogarithmic factors.
- Sequential communication: 2mn words make sequential TSQR optimal and a factor mn/(4W) fewer than blocked Householder QR.Sequential TSQR requires at most about 3mn/W messages, which is close to optimal.
- Sequential communication: 12mn^2/W^3/2 messages make sequential CAQR optimal, while the LAPACK and ScaLAPACK analyses assume W ≥ 2m.If that memory condition fails, those existing algorithms may communicate even more.
2. Tall-Skinny QR - TSQR.
TSQR computes the R factor through a reduction tree, preserving numerical stability while reducing communication for tall-skinny QR. Its parallel and sequential variants attain optimal communication costs and outperform stable alternatives asymptotically.
- TSQR as a reduction operation: TSQR repeatedly replaces stacked matrices with the R factor of their QR decomposition until one R factor remains.The combination pattern is a reduction tree, and any tree produces a valid QR decomposition.
- TSQR as a reduction operation: Any reduction tree yields a valid QR decomposition because the mapping from stacked matrices to their R factor is associative modulo roundoff and diagonal signs.The resulting Q factors can remain represented as a tree rather than being explicitly multiplied.
- Sequential TSQR: Sequential TSQR processes row blocks that fit in fast memory, retaining the accumulated R factor while moving Q representations between fast and slow memory.A flat tree supports this out-of-core procedure when fast memory holds slightly more than a fraction m/p of the rows.
- Parallel TSQR: Parallel TSQR performs log P messages and transfers 1/2 n^2 log P words along a binary-tree critical path.The critical path contains one local QR of an m/P × n block and log P QR factorizations of 2n × n stacked triangular factors.
- Comparison with alternatives: Parallel TSQR sends 2n times fewer messages than PDGEQRF while using the same bandwidth and remaining about as fast as Cholesky QR.The paper characterizes it as simultaneously fast and stable, unlike the compared alternatives.
- Comparison with alternatives: Sequential TSQR sends about mn/4W times fewer words and n/4 times fewer messages than PFDGEQRF, with fewer words than CholeskyQR when W ≥ n^2.TSQR is numerically stable because its computation consists of repeated orthogonal transformations.
3. Communication-Avoiding QR - CAQR.
CAQR combines right-looking QR with TSQR panel factorizations, choosing block and processor-grid parameters to reduce communication while controlling flop cost. Parallel CAQR matches ScaLAPACK in leading flops and words transferred but can require substantially fewer messages; sequential CAQR minimizes flops, bandwidth, and latency across fast-memory sizes.
- CAQR uses TSQR to factor panels within a right-looking QR algorithm, with the remaining work consisting of updates and bookkeeping.
- Parallel CAQR slightly increases flops only in lower-order terms while sending nearly the same number of words as ScaLAPACK’s PDGEQRF.
- A larger block size b can reduce ScaLAPACK’s 3n log Pr message term by a factor of b, but increasing b arbitrarily raises the flop count.
- With independently optimized b, Pr, and Pc, parallel CAQR matches ScaLAPACK in flops and words transferred while sending Θ(√(mn/P)) fewer messages.The comparison assumes n ≥ 2^11m/P, so the matrix is not too tall and skinny.
- Parallel CAQR’s bandwidth and latency terms match the corresponding lower bounds within constant factors under the stated parameter choices.
- Sequential CAQR minimizes flops, bandwidth, and latency for all fast-memory sizes W, unlike DGEQRF, RGEQRF, and RGEQR3.RGEQR3 attains the bandwidth lower bound but performs more floating-point operations than sequential Householder QR.
4. Lower Bounds for TSQR.
The section derives communication lower bounds for QR-related computations by analyzing smooth distributed functions and applies them to reduction-based matrix factorizations. These results support bandwidth lower bounds for CAQR and account for special handling of complex QR data.
- TSQR lower bounds: log P messages are both necessary and sufficient for parallel TSQR, achieved through a binary reduction tree.The lower bound follows from computing a nontrivial function of data distributed across P processors.
- Communication lower-bound framework: Smoothness restrictions make communication lower bounds meaningful by excluding lossless space-filling-curve compression of floating-point data.The analysis uses continuously differentiable functions and Jacobian conditions to characterize when fewer than all arguments can be communicated.
- Communication lower-bound framework: For a bijective function of P2’s n arguments, communicating fewer than n words is impossible; sending all n words is sufficient.The necessity follows because otherwise the function would be constant along a nontrivial set of P2’s inputs.
- Matrix-factorization lower bounds: c(c + 1)/2 words are necessary and sufficient for P1 to compute the Cholesky factor associated with P2’s contribution to the QR reduction.P2 can communicate either the upper-triangular entries of A2^T A2 or the entries of its Cholesky factor R.
- Matrix-factorization lower bounds: Complex QR requires adjusted counting because diagonal entries of R are generally real and each row has a unit-modulus complex scaling freedom.
5. Lower Bounds for CAQR.
The paper extends matrix-multiplication communication lower bounds to QR and LU, including latency, and proves that CAQR attains the resulting bounds within a constant factor.
- Lower-bound framework: Matrix-multiplication bandwidth lower bounds are extended to latency lower bounds and shown to apply to dense LU and QR decompositions.The latency bound follows because at most W words can move in one message, so a bandwidth lower bound B implies at least B/W messages.
- Lower-bound framework: The lower-bound analysis assumes conventional Θ(n3) matrix multiplication with commutative and associative reorderings, excluding asymptotically faster algorithms such as Strassen.Extending communication lower bounds to those faster algorithms remains open.
- Sequential CAQR: Sequential CAQR attains its communication lower bounds to within a constant factor.The result concerns both words moved between fast and slow memory and the number of messages.
- Parallel CAQR: Parallel CAQR attains its communication lower bounds to within a constant factor.The bounds apply to words sent and received and to messages involving at least one processor.
- Parallel CAQR: At least one processor must communicate a lower-bounded amount in parallel QR, with the square case covered when P ≥ 2^11.The square-case statement is given for the condition P ≥ 2^11.
- Arithmetic lower bounds: For columnwise QR, computing column j + 1 of R requires at least mj − j(j + 1)/2 multiplications, with additional bounds for later columns.The lower bound is based on the number of parameters needed to represent orthogonal matrices, whose Stiefel-manifold dimension is mj − j(j + 1)/2.
6. Related work.
Related work includes sequential and parallel variants of TSQR and CAQR, theoretical parallelism bounds, and prior matrix-multiplication communication bounds. The paper identifies parallel CAQR as novel and leaves sequential LU latency optimality open.
- TSQR and tiled QR: Prior work developed variations of sequential TSQR, often called tiled QR, and used them as panel factorizations or in shared-memory task-parallel QR.Those works were less focused on minimizing latency because communication costs were low in single-socket multicore settings.
- Parallel CAQR: Parallel CAQR is described as novel to the authors’ knowledge.Existing theoretical work instead studied exploitable parallelism and critical paths for QR based on Givens or Householder reductions.
- Communication lower bounds: Hong and Kung, followed by Irony, Toledo, and Tiskin, established communication lower bounds for sequential and parallel matrix multiplication.This paper attempts to extend those bounds to LU and QR factorization.
- LU factorization: The paper notes recent LU algorithms attaining at least some lower bounds, while sequential LU latency optimality remains an open question.One cited sequential LU result attains the bandwidth lower bound, but whether it attains the latency bound is unresolved.
7. Conclusions and Open Problems.
The paper argues that communication-optimal QR techniques may extend to hierarchical machines and broader linear-algebra routines, but several architectural and algorithmic obstacles remain open.
- Communication lower bounds may extend to other Θ(n3) dense linear-algebra routines, including eigenvalue problems.
- The proposed hierarchical extension targets machines with multiple memory and parallelism levels, each having different bandwidths and latencies.
- CAQR cannot use simple recursion alone because its subtasks mix smaller QR decompositions with matrix multiplication.
- Recursive decomposition must handle structured small QR problems, such as pairs of triangles, without incurring dominating overhead.
- The analysis assumes homogeneous machines, so asymmetric communication costs and specialized processors remain outside its scope.
- An adaptive autotuning approach is identified as likely necessary for these more complicated architectures, leaving the issue as future work.